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inside  the  reflector  and  D.F.=5.0  .... 
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The  flow  diagram  of  OUTERS  subroutine 
The  flow  diagram  of  DT  subroutine  . . 
The  flow  diagram  of  INNER 


remaining  portion 


diffusion  theory  is  then  used  for  the 
the  reflector. 

For  GCRs  with  uniform  fuel  gas  density  distributions  in 
the  core  and  a 100  cm  thick  BeO  reflector  the  hybrid  model 
gives  accurate  keff  values  and  flux  profiles  as  compared  to 
those  from  the  standard  sn  approximation.  For  4 energy 
groups,  the  hybrid  model  calculation  time  is  roughly  five 
times  faster  than  for  a standard  S„  calculation:  for  21 
energy  groups,  the  hybrid  model  is  about  five  and  a half 
times  faster  than  the  standard  Sn  approximation.  In 
general,  it  is  found  that  the  hybrid  model  point  of 
transition  from  s„  theory  to  diffusion  theory  can  be  located 
closer  to  the  core-reflector  interface  when  the  core  density 

Calculations  performed  on  configurations  with 
nonuniform  fuel  gas  density  distributions  in  the  core  show 
that  the  hybrid  model  gives  the  same  kt((  and  flux  profile 
results  as  the  standard  sn  approximation.  In  contrast, 
diffusion  theory  is  incapable  of  treating  this  situation  in 
the  gas  core.  For  these  configurations,  it  is  also  found 
that  the  hybrid  model  is  more  than  five  times  faster  than 
the  standard  Sn  approximation. 

When  the  reflector  thickness  is  increased  from  100  cm 
to  150  cm,  the  computation  time  advantage  of  the  hybrid 
model  compared  to  the  standard  Sn  approximation  is  increased 
from  a factor  of  five  to  seven.  For  a reflector  thickness 


200  cm,  the  hybrid  model  is  over  eight  times 
standard  Sn  approximation. 


INTRODUCTION 

Gas  Core  Reactors  (GCRs)  have  undergone  both 
theoretical  and  experimental  studies  at  the  University  of 
Florida.  A significant  number  of  neutronics  and  energetics 
analyses  have  been  conducted.  The  results  of  this  research 
indicate  that  GCR  concepts  have  some  distinct  advantages  for 
power  generation  compared  with  the  other  nuclear  power 
systems.  For  space  power  systems,  these  advantages  include 

A.  high  working  fluid  temperatures  which  are 
advantageous  for  power  cycle  efficiency  and 
especially  for  heat  rejection  in  space; 

B.  the  temperature  of  the  coolant/working  fluid 
(fuel)  in  the  core  can  be  considerably  hotter  than 
for  the  surrounding  structure; 

C.  the  simple  geometry/ structure  of  the  core  helps  to 
minimize  thermal  stress  and  thermal  shock; 

D.  improved  safety  on  launch  and  orbit  as  compared  to 
solid  fueled  systems; 

E.  rapid  start-up  capability;  and 

F.  pulsed  power  capabilities  which  can  reduce  power 
conditioning  requirements. 


disadvantages . 


corrosive  properties  of  the  fluoride-fuel  material  (UF6 
and/or  UF;)  and  erosion  effects  due  to  the  high  temperature 
and  high  velocity  flow  of  the  fuel  and  working  fluids,  some 
type  of  protective  coating  is  needed  on  the  wall  surfaces 

contact  with  the  fuel/working  fluid  (e.g.  heat  exchanger, 
nozzle,  and  MHD  generator) . Also,  a major  concern  is  the 
possible  condensation  of  uranium,  uranium  compounds,  or 
other  metallic  compounds  on  walls  where  the  wall 
temperature  is  below  the  condensation  temperature  (e.g.  the 
MHD  generator  walls) . 

summary  of  Conclusions  from  Previous  SCB 
Heutronic  studies  at  the  ,y Diversity  of  Florida 

A large  number  of  neutronic  studies  have  been  performed 
on  GCRs  at  the  University  of  Florida.  Some  of  these  studies 

flux  distribution  calculations.  These  studies  [1]  have 
found  standard  diffusion  theory  to  generally  be  invalid  in 
the  gas  core  region  since  the  gas  core  dimensions  rarely 
exceed  a mean  free  path.  Therefore,  in  general,  it  is 
necessary  to  use  higher  order  approximations  to  the 
transport  equation  (e.g.  Sn  transport  theory)  when  analyzing 

Analysis  of  a variety  of  concepts  led  to  some  general, 
but  not  infallible,  conclusions.  For  example,  for  kc|(  or 


critical  mass  calculations,  "for  spherical  or  square 
cylindrical  systems  which  possess  a significant  degree  of 
grayness  for  the  gas  core  region,  multigroup  diffusion 
theory  is  generally  adequate  provided  good,  fast  and  thermal 
constants  are  used  and  provided  there  is  a multithermal 
approach  which  allows  for  full  upscattering  and 
downscattering.  Because  of  the  wide  range  of  geometries, 
temperatures,  temperature  differences,  pressures,  densities, 
and  materials  that  can  be  employed  in  gas  core  reactors, 
each  new  concept  demands  individual  scrutiny.  Reactor 
physics  techniques  that  work  for  one  gas  core  concept  may 
not  be  adequate  for  another;  few  group  diffusion  theory 
calculations  suffice  for  some  systems,  while  multigroup  Sn 
transport  theory  calculations  are  required  for  others. "[1, 
p.  134]  "Although  for  spherical  or  square  cylindrical 
systems  which  posses  a significant  degree  of  grayness  one- 
dimensional, few  group  diffusion  theory  is  often  adequate 
for  kef)  and  critical  mass  calculations,  accurate 
determination  of  neutron  flux  distributions  generally 
requires  at  least  one-dimensional  Sn-  calculations. 

Reliable  spatial  distributions  are  to  be  obtained  from  few- 
group  calculations,  while  reliable  spectral  distributions 
will  require  multigroup  calculations." [1,  p.  135] 


In  this  study,  neutronics  calculations  are  performed  on 
the  GCR  in  one  dimension  to  determine  the  impact  of  the 
different  neutronics  models  on  the  predictions.  The  XSDRNPM 
[2)  code  solves  either  the  Sn  approximation  to  the  transport 
equation  or  the  diffusion  equation;  XSDRNPM  is  used  to 
perform  123  group  (30  thermal  groups)  calculations  on 
selected  GCR  configurations.  The  group  constants  are  then 
collapsed  to  26  group  (6  thermal  groups)  and  to  4 group  (1 
thermal  group)  structures  (see  Appendix  G for  more 
information  about  the  energy  group  structures  of  these 
different  cross-section  libraries) . These  four  group  cross- 
sections  are  used  to  perform  keff  calculations  and  obtain 
information  about  the  angular  and  scalar  neutron  flux 
behavior  for  different  UFt  densities  and  different  gas  core 
and  reflector  region  dimensions  for  the  homogeneous  core, 
externally-moderated  "cavity"  reactor  configuration  shown  in 

An  examination  of  the  results  of  these  studies 
indicates  that  the  neutron  flux  shows  a high  degree  of 
nonlinear  anisotropy  for  high  gas  core  densities;  a very 
weak  to  moderate  degree  of  nonlinear  anisotropy  for 
intermediate  densities;  and  very  weak  nonlinear  anisotropy 
to  essentially  linear  anisotropy  for  low  densities. 

However,  even  for  the  linearly  anisotropic  systems,  since 
the  gas  core  dimensions  rarely  exceed  a mean  free  path,  and 


because  the  £,/Et  ratio  is  very  small  (-0.2 
diffusion  theory  approximation  is  not  valid  within  the  gas 
core  region.  Therefore,  a higher  order  approximation 
(e.g.,  an  Sn  approximation)  should  be  used  in  the  gas  core 
region  of  the  GCR  systems. 

The  use  of  the  standard  Sn  approximation  for  the  GCR 


computation  time  as  compared  with,  for  example,  the  CPU  time 
of  the  diffusion  theory  approximation,  especially  in  the 
case  of  multigroup  and/or  multi-dimensional  and/or  thick 
moderator-reflector  problems.  Therefore,  the  primary  aim  of 
this  dissertation  is  to  investigate  and  develop  a 
computational  neutronics  model  for  the  generic,  externally 
moderated  type  of  GCR  illustrated  in  Figure  1-1  that 
maintains  good  accuracy  and  significantly  shortens  the  CPU 
computation  time  of  the  standard  Sn  approximation,  i.e.,  to 
develop  an  effective  acceleration  method  for  Sn  transport 
theory  calculations  on  GCRs.  While  CPU  times  required  for 
static  neutronic  calculations  can  be  long,  they  are  still 
acceptable;  what  is  not  acceptable  is  the  required  CPU  time 
for  space-time  dependent  dynamic  neutronic  calculations. 

The  developed  acceleration  method(s)  will  be  an  "enabling 
tool"  insofar  as  GCR  space-time  dependent  calculations  are 
concerned.  This  is  because  space  time  dynamic  calculations 
will  require  a prohibitively  long  CPU  time  if  currently 
available  acceleration  methods  in  Sn  transport  theory  codes 
are  employed  (e.g.,  the  estimated  CPU  time  required  to 


analyze  a 10-second  GCR  transient  with  the  CRAV  X-MP/48 
Supercomputer  is  more  than  30  hours  for  two-dimensional,  4- 
group  calculations,  and  more  than  1 hour  for  one 
dimensional,  4-group  calculations).  The  development  of 
effective  acceleration  methods  for  S„  transport  theory 
calculations  on  GCRs  will  significantly  reduce  the  required 


One  promising  approach  to  achieving  this  goal  is  to  use 
the  Sn  approximation  in  the  gas  core  region  and,  if 
necessary,  through  the  first  several  mean  free  paths  in  the 
reflector;  the  diffusion  theory  approximation  is  then  used 
in  the  rest  of  the  reflector  (i.e.  a "hybrid"  s„-dif fusion 
theory  model  is  employed) . The  use  of  the  diffusion  theory 
approximation  for  the  rest  of  the  reflector  can 
significantly  reduce  the  CPU  time  required  for  the  solution 
of  the  problem. 


The  performance  of  this  hybrid  model  may  then  be 
further  improved  by  using  vectorization  (vector  processing) 
and/or  multitasking  (parallel  processing) . Vector 
processing  is  a technique  for  reducing  the  amount  of 
processor  time  required  to  run  a program.  It  can  be  applied 
to  segments  of  a FORTRAN  code,  such  as  a Do-Loop,  that 


performs  the  same  sequence  of  operations  o 
elements  of  arrays.  Multitasking  is  designed  to  improve 
turnaround  time  for  certain  computer-intensive  applications. 
Multitasking  can  be  used  on  programs  that  can  be  divided 


independent 


multiprocessor  configuration.  Since  parallel  processing  has 
been  shown  to  speed  up  "standard"  S„  transport  theory 
calculations,  it  will  also  be  advantageous  to  employ 
parallel  processing  in  configuration  with  the  developed 
hybrid  model.  This  is  because  many  different  processors  can 
be  used  at  the  same  time,  one  for  each  energy  group.  The 
combination  of  these  two  acceleration  techniques  can  yield  a 
significant  reduction  in  the  CPU  computation  time.  If  the 
developed  hybrid  model  subroutine  is  written  and  arranged  in 
a way  that  takes  advantage  of  vector  processing,  then 
vectorization  could  also  be  advantageous  (see  Appendix  C for 
more  information  about  vectorization  and  multitasking) . 

The  acceleration,  in  this  work,  will  be  achieved 
primarily  by  developing  a computational  "hybrid"  neutronics 
model  for  the  generic,  externally-moderated  type  of  GCR  (see 
Figure  1-1)  that  maintains  good  accuracy  and  significantly 
shortens  the  CPU  computation  time  of  the  standard  Sn 
approximation . 


Brief  Description  of  the 

soUd-Cgre  Reactor  BenchaarK.Pmbj.sia 

A solid  core  reactor  with  an  H?0  reflector  is  studied 
as  a benchmark  problem  (see  Figure  1-2) . This  benchmark 
problem  was  selected  because  both  diffusion  and  Sn  transport 
theory  approximations  can  be  successfully  applied  to  this 
problem,  and  they  each  give  essentially  the  same  results. 
Hence,  a "hybrid"  model  (a  model  that  will  be  developed  in 


B.  SLAB  GEOMETRY 


this  work  to  obtain  reduced  computation  times  relative  to 
the  standard  S„  transport  theory  approximation  for  GCR 
calculations)  should  then  also  give  the  same  results  for 
this  system.  The  "hybrid"  model  is  to  use  the  S„ 
approximation  both  in  the  gas  core  region  and,  if  necessary, 
through  the  first  several  mean  free  paths  in  the  reflector; 
the  diffusion  theory  approximation  will  then  be  used  in  the 
rest  of  the  system.  Therefore,  this  benchmark  problem  is  to 
be  used  as  a diagnostic  tool  to  verify  whether  the  "hybrid" 
neutronics  model  that  will  be  used  in  the  system  of 
interest,  i.e.  the  gas  core  reactor,  has  been  properly 
implemented. 

For  this  solid  core  reactor  benchmark  problem  with  4 
energy  group  neutron  cross  sections,  the  diffusion  theory 
approximation  is  found  to  underpredict  keff  by  less  than  1% 
relative  to  the  S2Pfl  result,  while  it  overpredicts  kt(,  by 
less  than  1%  relative  to  the  S2P,,  StP3,  and  S6Ps  results  in 
slab  geometry.  These  calculations  have  been  repeated  for 
different  dimensions  and  geometries.  Comparison  of  neutron 
multiplication  factor  (keH)  values  indicates  that  for  this 
benchmark  problem,  diffusion  theory  provides  good  accuracy 
(<1%  error  in  kt,(} . An  examination  of  the  total  number  of 
iterations  and  the  CPU  computation  times  for  these 
calculations  shows  that  (as  expected)  diffusion  theory  is 
much  faster  than  the  s„  transport  theory  approximation. 


Coupled  neutron-gamma  calculations  on  GCR  shields  are 
also  performed.  The  purpose  of  the  shielding  studies  is  to 
obtain  preliminary  shielding  designs  for  a UFt-fueled  gas 

the  weight  and  volume  of  the  reactor  while  achieving  an 

In  this  study,  one-dimensional,  coupled  neutron-gamma  Sn 
transport  theory  calculations  are  performed  with  XSDRNPM  for 
a number  of  different  shielding  materials  and  dimensions. 
These  studies  are  conducted  in  order  to  find  shielding 
designs  that  meet  performance  specifications  for  neutron 
fluence  and  gamma  dose  to  protect  the  sensitive  spacecraft 
equipment.  Fluence  and  dose  scaling,  or  normalization,  for 
these  shielding  design  studies  is  achieved  by  assuming  a 
mission  time  of  one  hour  at  a core  power  level  of  1000  MWt 
(200  MWe) . Single-material  and  two-material  shielding 
analysis  has  been  performed  in  this  work.  The  hybrid  model 
developed  in  this  work  is  to  be  used  primarily  to  accelerate 
conventional  ke((  calculations  performed  using  sn  transport 
theory  codes  on  "standard"  GCR  configurations  (i.e., 
externally  moderated  cavity  reactors) . Due  to  the  tight 
coupling  between  the  core  and  reflector  regions  in  an 
externally  moderated  GCR,  the  gas  core/reflector  region 
calculations  can  be  significantly  affected  by  the  presence 
and  conversely,  the  shielding  neutronic 


shielding. 


behavior  is  strongly  affected  by  the  GCR  core/reflector 
combination.  Consequently,  shielding  calculations  for  the 
GCR  should  generally  be  performed  with  the  gas  core  and 
reflector  regions  present  (i.e.,  the  conventional 
“bootstrapping"  approach  to  the  performance  of  reactor 
shielding  calculations  should  generally  not  be  employed  for 
GCRs) . Hence,  the  developed  “hybrid"  acceleration  method 
which  is  to  be  used  for  conventional  keff  calculations  on 
"standard”  GCR  configurations  should  also  accelerate  GCR 
shielding  calculations  because  of  the  significant  reduction 
of  the  computation  time  needed  for  the  core/reflector  part 
of  the  problem. 


organi?atien-gt_the  pjssertatjpn 


Chapter  2 of  this  dissertation  discusses  the  highl 
of  previous  Gas  Core  Reactor  research  that  has  been 
performed.  This  chapter  includes  a historical  survey  a: 
brief  summary  of  the  methods  used  in  the  analyses  of  GC 
and  presents  key  results  obtained  from  these  previous 


Chapter  3 of  this  dissertation  discusses  the 
preliminary  calculations  performed  in  this  study  for 
neutronic  characterization  of  the  GCR  using  the  XSDRNPM 
computer  code.  These  studies  lead  to  a better  understanding 
of  the  neutronic  behavior  and  characteristics  of  the  GCR 


dissertation  includes  a detailed 


Chapter  4 of  this 
description  of  the  "hybrid"  S„-diffusion  theory  model 
developed  in  this  study  to  accelerate  Sn  GCR  calculations 
for  GCRs.  It  also  discusses  some  other  possible  approaches 
to  accelerate  Sn  calculations  for  GCRs. 

Chapter  5 includes  results  of  the  neutronic  studies 
performed  on  the  solid  core  reactor  benchmark  problem.  This 
chapter  also  includes  a detailed  description  of  the 
implementation  of  the  hybrid  model  into  XSDRNPM-S,  the 
SCALE4  version  of  XSDRNPM. 

Chapter  6 of  this  dissertation  presents  results  of 
neutronic  calculations  performed  on  the  GCR  using  the 
XSDRNPM-S  computer  code.  This  chapter  also  includes  and 
discusses  the  results  of  the  application  of  the  developed 
hybrid  model  to  the  GCR. 

chapter  7 covers  the  GCR  shielding  analysis.  These 
studies  are  conducted  on  a UF(-fueled  GCR  providing  burst 
mode  space  power  to  find  shielding  designs  that  meet 
performance  specifications  for  neutron  fluence  and  gamma 
dose  to  protect  the  sensitive  spacecraft  equipment. 

Chapter  8 draws  conclusions  from  these  results  and 
provides  some  recommendations  for  further  research  aimed  at 
improving  the  computational  efficiency  of  GCR  neutronics 
calculations. 

Appendix  A includes  a general  description  of  XSDRNPM-S, 
the  SCALE4  version  of  XSDRNPM. 


Appendix  B includes  the  subroutine  branching  structure 
for  the  XSDRNPM-S  code  and  a brief  description  of  the 
function  of  each  subroutine. 

Appendix  C includes  and  discusses  the  vectoriaation  and 
multitasking  methods  that  can  be  applied  to  improve  the 
computational  efficiency  of  the  developed  hybrid  model. 

Appendices  D,  E,  and  F include  general  descriptions  and 
flow  diagrams  of  the  OUTERS,  DT,  and  INNER  subroutines, 
respectively.  Subroutine  OUTERS  is  the  flow  control  routine 
for  the  flux  and  search  calculation.  It  calls  different 
subroutines  to  calculate  the  flux  according  to  the 
approximation  specified  by  the  user  and  it  loops  also  over 
energy  groups  that  define  an  outer  iteration  loop. 

Subroutine  DT  performs  the  diffusion  theory  flux  calculation 
required  on  the  first  outer  iteration  when  a fission  density 
guess  is  supplied,  or  when  requested  in  the  input  by  the 
user.  Subroutine  INNER  loops  over  the  angles,  calculating 
the  angular  flux  at  each  mesh  boundary  using  a weighted 
diamond  difference  scheme;  it  also  computes  updated  scalar 
flux  moments  during  the  sweep.  The  general  descriptions  and 
flow  diagrams  of  these  three  subroutines  help  to  provide  a 
better  understanding  of  how  the  developed  hybrid  model 
subroutine,  HYBRID,  is  used  to  provide  hybrid  model 
calculations  in  XSDRNPM-S. 

Appendix  G of  this  dissertation  contains  tables  which 
give  the  energy  group  structures  of  the  different  cross- 
section  libraries  used  in  this  work. 


CHAPTER 


GASEOUS  CORE  REACTORS-HISTORICAL  SURVEV 
PREVIOUS  STUDIES  AT  THE  UNIVERSITY  OF  FIORIDA 


Introduction 


n performed  on  many  different  ga 
e past  three  decades.  Several 
work  have  previously  been  per- 


Extensive  work  has 

literature  reviews  of  t 
formed. 

Gaseous  core  reactors  use  a gaseous  fissile  material  as 
fuel  for  power  generation.  The  fissioning  gas  "drives"  the 
reactor  core,-  a wide  variety  of  physical  arrangements  using 
external  or  internal  moderation  of  neutrons  and  methods  of 
containing  the  gas  and  extracting  the  energy  have  been 
suggested.  Scientific  research  has  been  conducted  on 
candidate  concepts  for  the  past  35  years.  Most  of  the 
research  has  been  limited  to  conceptual  phase  studies,  with 
a few  experimental  programs  conducted  to  determine  basic 
neutronics  and  fluid  characteristics  and  the  associated 
materials  and  materials  interaction  studies. 

of  nuclear  power  generation — including  improved  fuel 
utilization,  safety,  plant  efficiency,  special  high- 
performance  features,  load-following  capabilities  and  power 
conversion  optimization,  and  its  unique  potential  for  space 


power  applications — has  kept  the  gaseous  core  area  alive 
throughout  these  decades  of  erratic  directions  in  the 
development  of  nuclear  power. 

The  major  research  efforts  in  the  United  States  were 
coupled  to  the  nuclear  space  power  program;  these  were 
essentially  halted  in  1973  when  the  program  was  cancelled. 

during  the  1970s:  the  mixed-flow  gaseous  reactor  studies  of 
Los  Alamos  National  Laboratory  (LANL)  and  United  Technology 
Research  Center  (UTRC)  and  the  cyclic  gaseous  core  and 
heterogeneous  gaseous  core  reactor  studies  at  the  University 
of  Florida.  In  addition,  an  effort  to  determine  the 
potential  and  viability  of  gaseous  core  reactors  for  large 
commercial  electric  power  generation  was  conducted  under  the 
Nonproliferation  Alternative  Systems  Assessment  Program 
(NASAP) . 

Histgrical-S.gry.gx 


The  original  version  of  the  following  historical  survey 
belong  to  N.  J.  Diaz's  [3]. 

The  first  U.S.  analytical  studies  of  a gas  core  nuclear 
reactor  were  reported  by  Bell  [4]  of  LANL  in  1955.  The 
reactors  considered  consisted  of  spherical  cavities  fueled 
with  uranium  hexafluoride  (UFfi)  gas  and  surrounded  by  heavy 
water  (DjO),  beryllium,  or  graphite  moderating  reflectors. 
Other  early  analytical  studies  included  those  by  Safonov  [5] 
in  1958,  Ragsdale  and  Hyland  [6]  in  1961,  Hyland  et  al.  [7] 


I,  and  Herwig  and  Latham  (8) 


critical  particle  densities  of  fissionable  atoms  in  these 
systems  correspond  to  molecular  densities  of  gases  at  less 
than  atmospheric  pressure,  the  term  "cavity  reactor"  was 
applied  to  these  configurations. 

Configurations  investigated  included  U2S,- , U2K-,  and 
Pu2J, -fueled  systems;  moderating-reflector  materials 
considered  included  D20,  beryllium,  and  graphite.  The  early 
studies  were  directed  toward  determining  basic  reactor 
physics  parameters  for  these  systems.  Quantities  calculated 
included  critical  masses,  critical  sites,  material 
reactivity  worths,  temperature  reactivity  coefficients,  and 
potential  breeding  ratios.  The  results  included  evaluations 
of  the  relative  effectiveness  of  the  various  fuels  and 
moderating-reflector  materials. 

Analytical  studies  of  gaseous  core  nuclear  rocket 
engines  were  carried  out  by  Plunkett  [93  in  1967.  Latham 
[10-123  and  Latham  et  al.  [133  conducted  a series  of 
extensive  calculations  (1966-1969)  for  a nuclear  light-bulb 
(closed  system)  gaseous  core  reactor  rocket  engine.  Coaxial 
flow  (open  system)  gas  core  nuclear  rocket  engines  were 
examined  by  Hyland  [14)  in  1971.  The  light-bulb  and  coaxial 
flow  rocket  engine  systems  are  actually  plasma  core 
concepts,  since  the  fissioning  gases  are  assumed  to  be  fully 
ionized.  Parameters  calculated  again  included  critical 
masses,  critical  sizes,  material  worths,  pressure,  and 
temperature  reactivity  coefficients. 


comparison  ot  theoretical  predictions 


experimental  data  for  a gaseous  uranium  core  reactor  was 
performed  by  Kikoin  et  al.  [15]  and  reported  in  Moscow  in 
1959.  Neutronic  parameters  and  UP6  behavior  and 
conditioning  with  fluorinating  agents  were  reported.  The 
core  had  internal  moderation  (beryllium)  and  a graphite 
reflector.  It  went  critical  with  3340  g of  uranium  at  90% 
U23s . The  first  such  study  in  the  United  States  was 
conducted  by  Mills  [16]  in  1962.  In  1965,  Jarvis  and  Beyers 
[17]  of  LANL  made  a comparison  between  theoretical 
predictions  and  experimental  results  for  a D?0  reflected 
cavity  reactor. 

A series  of  critical  experiments  (1967-1969)  was 
performed  by  Pincock  et  al.  [16,19]  and  Kunze  et  al.  [20]  to 
test  the  ability  of  calculational  procedures  to  evaluate 
criticality  and  other  reactor  parameters  for  a configuration 
closely  resembling  a coaxial  flow  system.  A benchmark 
critical  experiment,  with  spherical  symmetry,  on  a gas  core 
nuclear  reactor  concept  was  reported  in  1972  by  Kunze  et  al. 
[21].  More  recent  experimental  research  programs  on  gas 


[22]. 


Much  of  the  early  gas  core  reactor  research  was 
motivated  by  a desire  by  the  U.S.  National  Aeronautics  and 
Space  Administration  (NASA)  to  apply  such  systems  for  rocket 
propulsion  in  space.  The  original  conceptual  gaseous  core 
nuclear  rocket  engines  are  specialized  systems  in  which 


nuclear  energy  is  generated  and  converted  into  thrust  by  the 
expulsion  of  heated  gases  in  the  plasma  state.  The 
application  of  gas  core  reactors  to  large  central  station 
electric  power  generation  also  had  some  early 
investigations.  In  1969,  Gritton  and  Pinkel  [23]  of  the 
Rand  Corporation  reported  one  such  study  in  which  a 4000-MW 
(thermal)  spherical  gaseous  core  system  with  a 5-ft  radius 
was  analyzed  at  a pressure  of  11  atm.  The  central  cavity 
was  surrounded  by  a moderating-reflector  region  and  by  banks 
of  energy  conversion  devices;  the  core  was  in  a plasma 

A significant  series  of  theoretical  and  experimental 
investigations  on  gaseous  reactors  was  conducted  by  Williams 
and  co-workers  [24-31]  at  The  Georgia  Institute  of 
Technology  from  1968  through  1975.  The  concepts 
investigated  centered  on  plasma  cores,  breeder  reactor  power 
plants,  and  advanced  energy  conversion  systems,  with 
extensive  work  on  magnetohydrodynamic  power  systems. 

A novel  concept  in  gas  core  reactor  design  was 
introduced  in  1974  by  Diaz  and  Dugan  [32,33]  and  Diaz  et  al. 
[34].  The  Heterogeneous  Gas  Core  Reactor  (HGCR)  contains  an 
array  of  moderator  channels  arranged  in  the  gaseous  core 
region.  The  HGCR  thus  has  "internal"  moderation  and  the 
core  heterogeneity  leads  to  significant  thermal-hydraulic 
and  energetic  advantages  along  with  improved  technological 
feasibility  with  respect  to  other  gas  core  reactor  concepts. 
The  HGCR  was  proposed  for  central  station  electric  power 


generation  and  utilizes  UP4  gas  at  temperatures  in  the 
vicinity  of  1000K,  in  contrast  to  the  plasma  core  concepts 
that  were  analyzed  at  significantly  higher  operating 
temperatures. 

Other  studies  of  gas  core  reactors  applied  to 
terrestrial  power  plants  were  reported  by  Rust  and  Clement 
[35]  in  1976.  They  examined  externally  moderated  UF4 
breeder  reactor  plants  in  which  the  blanket  was  molten  salt; 
nonionized  gases  were  used  for  the  fuel. 

A study  by  Lowry  [36]  in  1977  examined  the  Mixed-Flow 
Gas  Core  Reactor  (MFGCR)  concept.  The  MFGCR  employs  a 
cluster  of  seven  cavities  externally  moderated  by  beryllium 
and  surrounded  by  a molten  salt  breeding  blanket.  The  fuel 
gas  mixture  is  passed  through  the  core  region  by  circulators 
to  remove  heat  from  the  reactor  and  transport  it  to  the 
intermediate  heat  exchangers.  The  multicavity  concept 
provides  better  moderator  distribution  for  the  thermal 
reactor  than  a single  reflected  cavity  and  is  necessary  for 
the  mixture  temperature  control  required  in  the  mixed  flow 
concept.  "Mixed  flow"  refers  to  gas  mixing  at  the  outlet  of 
each  cavity;  a hot,  slowly  moving  portion  of  the  gas  in  the 
central  part  of  the  cavity  (average  exit  temperature  of  1000 
K)  mixes  with  cooler  swirl  flow  of  fuel  gas  that  maintains 
contact  with  and  cools  the  cavity  wall  to  800  K.  The  vortex 
separation  of  the  fuel  gas  into  two  streams  in  each  cavity 
is  accomplished  by  tangential  injection  slots  in  each 
cavity.  The  UF4  fuel  mixture  enters  a distribution  manifold 


cavities , 


at  the  top  of  reactor  vessel,  flows  through  the 
and  is  collected  in  a discharge  manifold  at  the  bottom  of 
the  vessel.  The  mixture  leaves  the  vessel  and  flows  through 
a heat  exchanger  and  a circulator.  A second  helium  loop 
cools  the  moderator  and  blanket.  An  experimental  program  on 

predicted  behavior  of  the  concept  [37].  The  program,  under 
NASA  sponsorship,  was  conducted  by  LANL  and  UTRC  at  Ijos 
Alamos.  Several  critical  configurations  were  successfully 
achieved  and  a nonfissioning  simulation  of  vortex  flow 
conditions  was  established. 

A significant  departure  from  the  above  steady-state  gas 
core  systems  is  to  be  found  in  two  pulsed  gas  core  reactor 
concepts.  The  Pulsed  Nuclear  Piston  (PNP)  system,  also 
known  as  the  nuclear  piston  engine,  was  first  studied  by 
Kylstra  et  al.  [38].  Early  studies  were  extended  by  Dugan 
[39]  in  1976  and  by  Diaz  et  al.  [40]  in  1978.  The  PNP 
consists  of  a small  pulsed  gaseous  core  reactor  enclosed  by 
a moderating-reflector  cylinder  and  piston  assembly  that 
operates  on  a thermodynamic  cycle  similar  to  the  internal 
combustion  engine.  A sister  concept  is  the  Pulsed  Gas 
Generator  (PCG)  reactor  proposed  by  Diaz  et  al.  [41,42]  in 
1979.  The  gas  generator  is  similar  to  the  nuclear  piston 
concept,  except  that  it  employs  a core  of  fixed  dimensions. 
The  PGG  concept  has  the  advantage  of  mechanical  simplicity, 
since  it  requires  no  moving  piston;  however,  the  PGG  concept 
has  a lower  efficiency  than  the  PNP  system  since  it  yields 


only  thermal  and  no  mechanical  or  shaft  power.  These 
systems  have  the  advantage  of  attaining  high  peak  values  of 
variables,  such  as  temperature  and  pressure,  while 
maintaining  relatively  low  cycle-averaged  values.  Power 
generation  for  these  pulsed  systems  per  chamber  is  of  the 
order  of  a few  megawatts  for  steady-state  use  and  up  to 
hundreds  of  megawatts  for  transient  operation.  Analytical 
studies  on  these  pulsed  cores  were  compared  with 
experimental  results  obtained  at  the  Plasma  Core  Assembly  at 
LANL  [42]. 

Extensive  analytical  studies  have  been  performed  on 
many  different  gas  core  reactor  concepts  over  the  past  30 
years.  Neutronic  experiments  conducted  during  this  period 
have  been  mostly  restricted  to  verifying  predictions  for 
criticality  and  basic  reactor  physics  parameters,  such  as 
reactivity  coefficients  and  neutron  lifetimes;  space-, 
energy-,  and  time-dependent  experiments  were  completed  at 
LANL  in  1980  [39].  Basic  materials,  fluid  flow,  and  heat 
transfer  experiments  have  also  been  conducted  for  some  of 
the  gas  core  reactor  concepts.  Gas  core  fission  concepts 
have  a series  of  attractive  features,  including  improved 
resource  utilization,  adaptability  to  different  energy 
conversion  systems,  and  good  control  and  safety 
characteristics  that  warrant  their  continued  development. 
Pulsed  or  cyclic  gaseous  core  nuclear  systems  have  the 
additional  advantages  of  compactness,  high-power  density, 
and  high  efficiency  and  appear  to  offer  a viable  energy 


range  of  applications. 


generation  alternative  for  a wide 
including  space-  and  land-based  power  units.  The 
establishment  of  a sound  scientific  data  base  that  includes 
key  information  in  the  area  of  neutronios,  thermophysical 
properties,  and  materials  for  cyclic  gaseous  core  reactors 
has  been  the  objective  of  a lengthy  theoretical/experimental 
research  program  at  the  University  of  Florida  [39-42]. 

Previous  Gas  Core  Reactor  Studies  at  the 


At  the  University  of  Florida,  research  on  Gaseous/Vapor 
Core  Reactors  (GCR/VCR)  continued  over  the  past  two  decades. 
During  these  years  many  different  concepts  were  studied  and 
characterized  in  the  Department  of  Nuclear  Engineering 
Sciences  and  the  Innovative  Nuclear  Space  Power  and 
Propulsion  Institute  (INSPI)  [43-73].  The  major  concepts 
that  have  evolved  from  these  research  efforts  are 
a)  law  power  Pulsed  Gaseous  core  Resotor 

The  Low  Power  Pulsed  Gaseous  Core  Reactor 
(PGCR)  was  one  of  the  first  GCR  concepts  studied 
at  the  University  of  Florida  [39].  The  (PGCR) 
consists  of  a cylindrical  gaseous  core  enclosed  by 
a moderating-reflector.  The  gaseous  fuel/working 
fluid  consists  of  a mixture  of  highly  enriched 
uranium  hexafluoride  (UF6)  and  helium  gas;  the 
helium  is  added  to  enhance  the  thermodynamic, 
transport  and  heat  transfer  characteristics  of  the 


cyclically 


primary  working  fluid.  Gaseous  fuel  is 
injected  into  the  core,  the  core  is  pulsed  to 
supercriticality  and  then  the  gas  is  exhausted. 

The  high  temperature  UFj/He  exhaust  gas  is  cooled 
in  a UFj/He-to-He  heat  exchanger.  The  heated 
helium  can  then  be  passed  either  directly  through 
gas  turbine  or  be  used  in  a steam  generator  to 
drive  a turbine.  This  pulsed  system  has  the 
advantages  of  attaining  a high  peak  gas 
temperature  (-1700  K)  and  high  thermal  neutron 
fluxes  (-1015  n/cm2-sec)  while  maintaining  low 
average  temperature  (-800  K) . 
b)  HeteraaenegHS..Si§-Core  Reactor 

The  core  heterogeneity  in  the  Heterogeneous 
Gas  Core  Reactor  (HGCR)  improves  the  heat  transfer 
characteristics,  the  power  and  temperature 
distributions  and  neutron  and  fuel  economy 
compared  to  the  previous  gaseous  core  reactor 
systems  investigated  at  the  University  of  Florida 
[33] . The  HGCR  consists  of  a reactor  pressure 
vessel,  core  barrel,  a gaseous  core  with  its  array 
of  moderator  channels  or  cells  and  a reflector 
region.  The  gaseous  fuel  (a  mixture  of  UF4-He) 
flows  through  the  core  and  transports  a major 
portion  of  the  thermal  energy  generated  within  the 
core  by  the  fissioning  gaseous  fuel  to  an  outside 
heat  exchanger  where  its  heat  is  transferred  to  a 


secondary  coolant  which  passes  through  a turbine. 
An  auxiliary  coolant,  such  as  helium,  can  be  made 
to  circulate  through  the  channels  in  the  moderator 
to  provide  additional  cooling,  if  necessary.  HGCR 
systems  with  blanket  regions  surrounding  the  core 
have  also  been  examined. 

Neutronics  and  energetics  studies  of  the  HGCR 
indicate  that  it  offers  the  possibility  of 
significant  gains  in  overall  thermal  efficiency 
over  Heavy  Hater  Reactors  (HWR's)  and  Light  Water 
Reactors  (LWR's).  The  high  temperature  gas  flow  in 
the  HGCR  allows  the  reactor  to  interface  with 
proven  high  efficiency  power  conversion  systems. 
The  Bimodal  Gaseous  Core  Reactor  Power  System 

Investigations  have  also  been  performed  on 
the  Bimodal  Gaseous  Core  Reactor  (BGCR)  Power 
System  (68,69].  The  BGCR  consists  of  two  reactor 
systems,  namely,  the  central  high  (burst)  power 
gas  core  reactor  (BPGCR)  system  which  is 
surrounded  by  an  annular  ring  of  small, 
cylindrical  low  power  PGCR  chambers  embedded  in 
beryllium  moderator.  The  identical  upper  and 
lower  halves  of  the  reactor  system  are  separated 
by  a common  moderator  slab  at  the  mid-plane. 

The  annular  ring  of  PGCR  chambers  is  designed 
to  provide  low  power  for  station  keeping/ 
surveillance  purposes.  The  primary  working  fluid 


(fuel)  is  a mixture  of  uranium  hexafluoride  (uf6) 
and  helium  gas.  During  low  power  operation  of  the 
system,  the  central  BPGCR  chamber  will  either  be 
voided  or  have  a constant,  low  pressure  nonfuel 

The  transition  from  station  keeping  to  the 
high  or  burst  power  mode  of  operation  can  be 
achieved  by  circulating  a gaseous  uranium  fuel 
through  the  central  BPGCR  cavity.  The  neutronic 
coupling  effects  of  the  PGCR  chambers  on  the 
gaseous  fuel  bearing  central  chambers  provide  the 
necessary  criticality/heating  condition  to 
generate  a partially  ionized  plasma  {=  2000-4000 
K) . The  high  pressure  fuel  gas  mixture  is  driven 
into  a disk  MHD  generator  through  a supersonic 
nozzle.  The  energy  conversion  occurs  in  the  disk 
MHD  generator.  The  fuel  gas  then  passes  through  a 
diffuser  and  a radiator  heat  exchanger  before 
being  circulated  back  into  the  central  BPGCR 
chamber  by  a compressor. 

In  addition  to  providing  station  keeping 
power,  the  PGCRs  relax  the  reactivity  requirements 
for  the  central  BPGCR.  This  provides  greater 
flexibility  in  selecting  the  BPGCR/MHD  system 
fuel/working  fluid  composition.  Since  the  PGCRs 
provide  a base  level  neutron  flux  in  the  system 
and  maintain  the  system  "hot",  a more  rapid 


transition  to  burst  power  is  possible  than  if  the 

Ultrahioh  Temperature  Vapor  Core  Reactor 

The  current  focus  of  GCR/VCR  research  efforts 
at  INSPI  is  the  burst  power  Ultrahigh  Temperature 
Vapor  Core  Reactor  (UTVR)  using  MHD  energy 
conversion  [64,65].  The  UTVR/MHD  system  employs 
uranium  tetrafluoride  (UFt)  as  the  fissile  fuel 
and  an  alkali  metal  fluoride  (KF,  NaF  or  LiF7)  as 
the  working  fluid. 

The  UTVR  consists  of  a central,  cylindrical 
core  which  contains  a vapor  mixture  of  UF4  fuel 
and  the  metal  fluoride  and  is  surrounded  by  a few 
boiler  columns  in  which  the  UF4  liquid  is 
vaporized  before  entering  the  central  UTVR  core 
region.  The  central  core  region  and  the  boiler 
regions  are  embedded  in  a beryllium  oxide 
reflector/moderator.  After  exiting  the  core,  the 
vapor  fuel/working  fluid  mixture  passes  through  a 
nozzle,  disk  MHD  generator  and  diffuser.  The 
fuel/working  fluid  mixture  is  then  condensed  in  a 
condensing  radiator/heat  exchanger.  The  condensed 
UF4  is  pumped  back  into  the  boiler  while  the 
condensed  metal  fluoride  serves  as  a wall  coolant 


CHAPTER 


theory  equation;  XSDRNPM  is  used  to  perform  123  group  (30 
thermal  groups)  calculations  on  selected  GCR  configurations. 
The  group  constants  are  then  collapsed  to  26  group  (6 
thermal  groups)  and  to  4 group  (1  thermal  group)  structures. 
These  4 group  cross-sections  are  used  to  perform  kBff 
calculations  and  obtain  details  about  the  angular  and  scalar 
neutron  flux  behavior  in  the  GCR  for  different  UF4  densities 
and  different  gas  core  and  reflector  dimensions. 

For  the  spherical  gas  core  reactor  system  shown  in 
Figure  1-1,  examination  of  the  predicted  fast  and  thermal 
group  flux  distributions  shows  several  important  trends. 

The  thermal  and  fast  neutron  flux  distributions  versus 
radial  position  are  shown  in  Figures  3-1  and  3-2, 
respectively,  as  a function  of  fuel  atom  density.  In  the 
gas  core  region,  the  thermal  neutron  flux  level  decreases 
with  increasing  UFt  density  due  to  the  increasing 
absorption.  The  fast  neutron  flux  in  the  gas  core  region 
increases  with  increasing  gas  loading  due  to  increased 
neutron  production  from  fission. 

As  also  shown  in  Figure  3-1,  the  thermal  neutron  flux 
peak  in  the  Beo  reflector/moderator  increases 
with  increasing  fuel  atom  density  in  the  core.  This 
increase  is  due  to  the  greater  relative  availability  of  fast 
neutrons  in  the  core  which  then  slow  down  in  the  reflector. 
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The  increased  fast  neutrons  in  the  core  are  due  to  the 
increased  fuel  atom  density  and  fission  rate. 

Comparative  neutronics  calculations  are  performed  on 
the  GCR  in  one  dimension  to  determine  the  impact  of 
different  neutronics  models  on  the  predictions.  Various  UFt 
gas  loadings  and  various  reflector  thicknesses  are  used  as 

total  and  scattering  mean  free  paths  with  gas  loading.  This 
table  also  shows  that  the  dimension  of  the  gas  core  exceeds 
the  mean  free  path  of  the  neutrons  only  for  very  high 
density  factors.  Finally,  from  this  table,  the  ratio  of 
s,/st  ^ or  thermal  neutrons  in  the  gas  core  is  found  to  be 
very  small  (-0.2).  These  considerations  indicate  that  the 
diffusion  theory  approximation  within  the  gas  core  region  is 

The  kt11  results  from  StP5  and  diffusion  theory 
calculations  performed  with  XSDRNPM  are  compared  with  one 
another  in  Table  3-2.  This  table  shows  that  the  diffusion 
theory  approximation  is  found  to  overpredict  kg(f  by  less 
than  1%  relative  to  an  StP3  approximation.  An  examination 
of  the  total  number  of  iterations  and  the  calculation  times 
in  this  table  shows  that  diffusion  theory  is  much  faster 
than  the  Sn  transport  theory  approximation  (a  factor  of  up 
to  20  times  faster  for  4 groups  and  a factor  of  up  to  16 
times  faster  for  21  groups) . The  results  in  Table  3-2 
indicate  that  such  relative  savings  in  computation  time  are 
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more  significant  for  thick  reflectors. 

The  diffusion  theory  calculations  for  the  gas  core 
reactor  were  found  to  converge  in  just  a few  iterations 
(e.g.,  only  two,  or  at  most,  three  iterations  are  generally 
needed  to  reach  essentially  total  convergence) . In 
contrast,  the  S approximation  calculations  generally 
required  a significant  total  number  of  iterations  (e.g., 
greater  than  40)  to  reach  a convergence  level  of  about  10  . 
For  both  the  diffusion  theory  and  S„  calculations,  a flat 
fission  density  guess  was  used  to  start  the  calculations. 
Because  the  diffusion  theory  calculation  converges  very 
rapidly  in  this  case,  comparison  of  computation  time  with 
that  of  the  corresponding  Sn  calculation  can  be  somewhat 
misleading.  Therefore,  a bad  initial  flux  guess  is  used  for 
the  diffusion  calculations  in  order  to  increase  the  number 
of  iterations  and  calculation  time.  In  order  to  fairly 
compare  the  S„  results  with  the  diffusion  theory  results, 
the  same  bad  guess  is  then  also  used  for  the  corresponding 
Sn  calculation. 

The  angular  thermal  neutron  flux  behavior  at  the 
interface  between  the  core  and  the  reflector  and  at  the 

reflector  interface,  as  a function  of  fuel  atom  density,  are 
shown  in  Figures  3-3  to  3-12.  An  examination  of  these 
figures  indicates  that  the  neutron  flux  exhibits  an 
essentially  linear  anisotropy  to  a very  weak  nonlinear 
anisotropy  for  D.F.=0.5,  D.F.=1.0,  and  D.F.«3.0;  a weak  to 
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Thermal  angular  neutron  flux  behavior  of 
the  spherical  GCR  at  the  interface  between 
the  core  and  the  reflector  (D.F.=3.0). 


D.F=5.0 


Thermal  angular  neutron  flux  behavior  of 
the  spherical  GCR  at  the  interface  between 
the  core  and  the  reflector  (D.F.S5.0). 
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- The  flux  has  been  normalized  to  1 (fission/s) 

- D.F. =100.0  corresponds  to  the  fuel  atom  densities 


Thermal  angular  neutron  flux  behavior  of 
the  spherical  GCR  at  the  interface  between 
the  core  and  the  reflector  (D.F. =100.0) . 


moderate  degree  of  nonlinear  anisotropy  for  D.F.=5.0  and 
D.F. =8.0;  and  a high  degree  of  nonlinear  anisotropy  for 
D.F. =80.0  and  D.F. =100.0.  However,  even  for  the  linearly 
anisotropic  systems,  since  the  gas  core  dimensions  rarely 
exceed  a mean  free  path,  and  because  the  £,/£(  ratio  is  very 
small  (-0.2),  the  diffusion  theory  approximation  within  the 
gas  core  region  is  not  valid.  Therefore,  a higher  order 
approximation  (e.g.,  an  Sn  approximation)  should  be  used  in 
the  core  region  of  the  GCR.  On  the  other  hand,  Figures  3-13 
and  3-14  show  the  angular  thermal  neutron  flux  behavior  at 
different  points  inside  the  reflector  when  D.F. =100.0  in  the 
core.  These  results  indicate  that  the  neutron  flux  will  be 
nonlinearly  anisotropic  up  to  a certain  distance  into  the 
reflector  (Beo  reflector) , which  is  found  to  be  about  5 to  6 
cm;  after  this  point  a diffusion  theory  approach  can  be  used 
instead  of  a higher  order  approximation.  Of  course,  one 
would  expect  that  this  distance  into  the  reflector  will  be 
smaller  for  density  factors  lower  than  100.0.  This 
particular  density  factor  (i.e.  D.F. =100.0)  has  been  used  in 
Figures  3-13  and  3-14  because  it  has  the  highest  degree  of 
nonlinear  anisotropy  of  the  examined  systems.  These  results 
indicate  that  it  should  not  be  necessary  to  extend  the  S„ 
approximation  beyond  the  first  5 to  6 cm  of  the  BeO 
reflector.  The  use  of  a diffusion  theory  approach,  instead 
of  a higher  order  approximation  (e.g.,  an  s„  approximation, 
which  can  require  long  calculation  times  as  compared  with 


- The  flux  has  been  normalized  to  1 (fission/s) . 

- 1.  for  a point  at  the  interface  between  the  core 

and  the  reflector. 

- 2.  for  a point  at  3.33  cm  inside  the  reflector. 

- 3.  for  a point  at  6.66  cm  inside  the  reflector. 

- 4.  for  a point  at  9.99  cm  inside  the  reflector. 

- 5.  for  a point  at  13.33  cm  inside  the  reflector. 


Figure  3-13 


Relative  angular  thermal  neutron  flux 
behavior  at  different  points  inside  the 
reflector  (D.F.-100. 0) . 
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the  diffusion  theory  approximation)  in  the  rest  of  the 
reflector  can  significantly  reduce  the  computation  time,  as 
can  be  inferred  from  the  results  in  Table  3-2. 

In  addition,  for  the  spherical  gas  core  reactor  system 
shown  in  Figure  1-1,  examination  of  the  predicted  thermal 
neutron  flux  behavior  from  both  the  Sn  and  diffusion  theory 
approximations  shows  some  more  important  features.  Figure 
3-15  shows  that  both  the  s„  and  diffusion  theory 
approximation  yield  almost  the  same  thermal  neutron  flux 
behavior  throughout  the  reflector  region  while  there  is  some 
difference  within  the  gas  core.  However,  actual  gas  core 
reactors  will  have  a nonuniform  gas  density  distribution  due 
to  the  nonuniform  temperature  distribution  inside  the  gas 
core  (the  temperature  close  to  the  center  is  higher  than  the 
temperature  in  a region  close  to  the  reflector) . Therefore, 
gas  core  reactors  with  two  different  nonuniform  gas  density 
distributions  are  examined  to  determine  the  impact  of  the 
density  variation  on  the  flux  behavior  as  predicted  by  both 
Sn  and  diffusion  theory  approximations.  In  the  first  case, 
a constant  density  factor  of  1.0  is  used  in  the  central 
region  of  the  core  and  a constant  density  factor  of  10.0  is 
used  in  the  outer  region  of  the  core.  Figure  3-16  shows  the 
thermal  neutron  flux  distributions  from  both  diffusion 
theory  and  the  Sn  approximations  for  this  nonuniform  gas 
density  distribution.  This  figure  shows  that  the  sn  flux 
distribution  prediction  is  changed  relative  to  that  for  the 
corresponding  system  with  a uniform  gas  density  distribution 


Radial  Position  (cm) 


- The  flux  has  been  normalized  to  1 (fission/s). 

- D.F.=10.0  correspond  to  the  fuel  atom  densities 
given  in  Table  3-1. 

Figure  3-15  Thermal  neutron  flux  v ariation  versus 

radial  position  for  spherical  GCR  for  both 
diffusion  theory  and  Sn  approximation 
calculations  (D.F.=10.0). 
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- D. F.=l . 0 and  10.0  correspond  to  the  fuel  atom 
densities  given  in  Table  3-1. 


Figure  3-16.  Thermal  neutron  flux  variation  versus 
radial  position  for  spherical  GCR  for 
both  diffusion  theory  and  Sn  approxi- 
mation calculations  using  two  different 
density  factors  in  the  core  according 


throughout  the  core.  In  contrast,  the  flux  distribution 
predicted  by  diffusion  theory  remains  unaffected  by  the 
nonuniform  gas  density  distribution  (i.e.,  it  is  the  same  as 
for  a uniform  gas  density  distribution  throughout  the  core) . 
This  is  another  example  of  the  limitations  of  diffusion 
theory  when  treating  the  gas  core.  Thus,  it  is  apparent 
that  reliable  prediction  of  the  flux  distribution  in  the  gas 
core  requires  S„  transport  theory  calculations  or  some 
comparable  high  level  approximation.  Reliable  prediction  of 
the  flux  spatial  distribution  is  crucial  for  space  time 
dynamic  neutronic  calculations  and  for  eventually  being  able 
to  determine  the  degree  of  coupling  which  exists  between  the 
gas  density  field  and  the  neutron  field  in  a gas  core 

For  the  second  nonuniform  gas  density  distribution 
case,  a density  factor  of  10.0  is  used  in  the  central  region 
of  the  core  and  a density  factor  of  1.0  in  the  outer  region 
of  the  core.  Figure  3-17  compares  the  thermal  neutron  flux 
variation  versus  radial  position  from  both  S„  and  diffusion 
theory  approximations  for  the  indicated  nonuniform  density 
distribution  in  the  core.  This  figure  also  shows  that 
diffusion  theory  is  completely  incapable  of  treating  this 
situation  in  the  gas  core. 
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Radial  Position  (cm) 


- The  flux  has  been  normalized  to  1 (fission/sl. 

- D.F.=10.0  and  1.0  correspond  to  the  fuel  atom 
densities  given  in  Table  3-1. 


Figure  3-17 


Thermal  neutron  flux  variation  versus 
radial  position  for  spherical  6CR  for 
both  diffusion  theory  and  Sn  approxi- 
mation calculations  using  two  different 
density  factors  in  the  core  according 


CHAPTER 


DESCRIPTION  OF  THE  "HYBRID"  Sn-DIFFUSION  THEORY  MODEL  AND 
SOME  OTHER  POSSIBLE  APPROACHES  TO  ACCELERATE 
Sn  CALCULATIONS  FOR  GCRs 


identified 


Hybrid  Model 


one  very  promising  approach  is  to  use  the  Sn 
approximation  in  the  gas  core  and,  if  necessary,  through  the 
first  several  mean  free  paths  in  the  reflector.  The 
diffusion  approximation  will  then  be  used  in  the  rest  of  the 
reflector  (i.e.,  a "hybrid"  Sn-diffusion  theory  model  is  to 
be  employed) . The  use  of  the  diffusion  theory  approximation 
for  the  rest  of  the  reflector  will  reduce  the  CPU 
computational  time.  A great  deal  of  effort  can  be  wasted  if 
unnecessary  Sn  calculations  are  used  in  the  rest  of  the 
reflector. 

There  are  many  computer  codes  which  can  perform  either 
Sn  or  diffusion  theory  approximation  calculations  (e.g. , 
XSDRNPM-S,  xsdrnpm,  onedant,  ONETRAN,  etc..).  However,  none 
of  these  codes  has  the  ability  to  selectively  implement 
different  approximations  in  different  regions  of  the  system 
or  different  approximations  in  one  problem  (i.e.,  either 
diffusion  theory  is  used  everywhere  or  Sn  theory  is  used 
everywhere  in  a given  iteration) . The  hybrid  model  requires 
the  use  of  the  two  different  approximations  in  one  problem. 
Therefore,  one  of  the  computer  codes  mentioned  above  will 
require  modification  to  achieve  this  combination  of  both  the 
s„  and  diffusion  approximations  in  one  problem  (XSDRNPM-S  is 
modified  for  this  purpose;  for  details  of  this  modification 
see  chapter  5) . In  the  Sn  approximation,  the  angular 
variable  is  placed  so  to  speak  on  equal  footing  with  the 


spatial  variable  and  is  treated  explicitly;  discretization 
of  the  angular  variable  is  accomplished  by  using  the 
discrete  ordinates  approach.  In  contrast,  when  the  flux  is 
calculated  under  the  diffusion  theory  approximation,  there 
is  no  explicit  treatment  or  discretization  of  the  angular 
variable;  it  is  assumed  that  the  flux  can  be  adequately 
represented  by  only  a linearly  anisotropic  angular 
dependence.  To  implement  the  hybrid  model,  it  is  necessary 
to  define  suitable  boundary  conditions  at  the  physical 
and/or  mathematical  interfaces  between  these  two 
approximations  in  order  to  correctly  preserve  continuity  of 
the  flux.  A number  of  approaches  can  be  used  to  reach  the 
goal  of  hybrid  model  implementation;  four  of  these 
approaches  are  explained  below. 

1.  The  first  approach  (the  approach  used  in  this 

dissertation)  to  implement  the  hybrid  model 
depends  on  the  use  of  an  initial  diffusion  theory 
calculation  for  both  the  gas  core  and  reflector. 
Then  the  S„  approximation  is  used  for  the  gas  core 
region  and,  if  necessary,  for  the  first  few  mean 
free  paths  inside  the  reflector.  For  this 
approach,  it  is  necessary  to  define,  in  addition 
to  the  physical  interface  at  point  IMP,  the 
mathematical  interface  which  is  shown  as  point 
IMMP  in  Figure  4-1.  The  mathematical  interface  is 
a point  inside  the  reflector  where  the  angular 
flux  behavior  is  no  worse  than  linearly 
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Second  step: 


converged  diffusion  calculation 


Albedo , obtained  from  3*  and  J' 

, calculation  from  1=1  to  I=IMMP. 

Albedo  boundary  condition  employed 


Last  step:  The  final  flux  solution  is  generated  as  shown  below 


Figure  4-1  The  first  two  steps  and  the  final  flux  result 
in  the  hybrid  model  implementation  according 
to  the  first  approach. 


anisotropic  (i.e.,  where  diffusion  theory  is 
valid) . The  steps  followed  in  this  approach  are 
explained  as  follows: 

a.  As  shown  in  Figure  4-1,  the  diffusion 

theory  approximation  is  applied  in  both  the 
core  and  the  reflector  regions  (i.e.,  from 
1=1  to  I=IM) . In  this  diffusion  calculation, 
standard  boundary  conditions  are  used  and  the 
converged  diffusion  theory  flux  solution  for 
each  energy  group  in  the  core  is  used  as  a 
flux  guess  to  initiate  the  following  Sn 
calculation  in  the  core:  the  diffusion  flux 
at  I=IMMP  (®°|WP)  is  used  to  generate  the 
partial  current  values,  J*  and  J*,  at  the 
mathematical  interface  (I=IMMP)  for  each 
energy  group.  These  partial  currents  are 
then  used  to  generate  the  values  of  the 
albedo,  es|mp,  at  the  mathematical  interface 
for  each  energy  group  "g"  which  are  then  used 
as  a boundary  condition  at  I=IMMP  for  the  sn 
calculation  on  the  next  step.  The  albedo  is 
defined  in  this  work  as  a’,*,,  = (J'/J’)’llw. 

The  XSDRNPM-S  computer  code  had  to  be 
modified  to  generate  the  partial  current 
values  from  the  flux  and  to  calculate  albedos 
from  the  partial  currents  at  I=IMMP.  The 
diffusion  calculation  process  in  the  first 


step  is  continued  until  the  standard 
diffusion  convergence  criteria  are  meet. 

In  this  step  an  S„  calculation  is 
performed  in  the  core  region  and,  if  needed, 
through  the  first  few  mean  free  path  of  the 
reflector  (i.e. , from  1=1  up  to  the 
mathematical  interface,  IMMP) . The  albedos 
generated  from  the  diffusion  theory  solution 
at  the  mathematical  interface  are  used 

as  a boundary  condition  for  this  Sn 
calculation.  This  Sn  calculation  process  is 
continued  until  the  standard  Sn  convergence 
criteria  are  meet.  The  converged  Sn  flux 
solution  from  1=1  to  I=IMMP  {tp Sn(1  IWPJ)  is 


The  final  flux  generated  from  the  hybrid 
model  calculation  is  composed  of  two  kinds  of 
solution  fluxes.  These  include  the  Sn 
approximation  flux  from  1=1  to  I=IMMP 
(^sn(1  ,Wj)  generated  in  step  b,  and  the 
diffusion  theory  flux  from  I=IMMP  to  I=IM 
(^,i„  , generated  in  step  a,  as  shown  in 

The  second  approach  to  implement  the 
hybrid  model  also  depends  on  the 
use  of  an  initial  diffusion  theory  calculation  for 


reflector.  Then 


approximation  is  used  for  the  gas  core  region  and, 
if  necessary,  for  the  first  few  mean  free  paths 
inside  the  reflector.  For  this  approach  it  is 
also  necessary  to  define,  in  addition  to  the 
physical  interface  at  point  IMP,  the  mathematical 
interface  which  is  shown  as  point  IMMP  in  Figure 
4-2 . The  mathematical  interface  is  a point  inside 
the  reflector  where  the  angular  flux  behavior  is 
no  worse  than  linearly  anisotropic  (i.e.,  where 
diffusion  theory  is  valid) . The  steps  used  in 
this  approach  are  explained  as  follows: 
a.  As  shown  in  Figure  4-2,  the  diffusion 

theory  approximation  is  applied  in  both  the 
core  and  the  reflector  regions  (i.e.,  from 
I«1  to  I-IM) . In  this  diffusion  calculation, 
standard  boundary  conditions  are  to  be  used 
and  the  converged  diffusion  theory  flux 
solution  in  the  core  is  used  as  a flux  guess 
to  initiate  the  Sn  calculation  in  the  next 
step;  the  diffusion  flux  at  X=IMMP  (dD]wt,)  is 
used  as  a boundary  condition  for  the  Sn 
calculation. 

The  diffusion  calculation  process  in 
this  first  step  is  continued  until  the 
standard  diffusion  convergence  criteria  are 
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Last  step:  The  final  flux  solution  is  generated  as  shown  below 


e second  two  steps  and  the  final  flux  result 
the  hybrid  model  implementation  according 
the  second  approach. 


In  this  step  an  Sn  calculation  is 
performed  in  the  core  region  and,  if  needed, 
through  the  first  few  mean  free  paths  of  the 
reflector  (i.e.,  from  1=1  up  to  the 
mathematical  interface,  IMMP) . The  diffusion 
theory  flux  solution  at  the  mathematical 
interface  (0°,w)  calculated  above  in  step  a 
is  used  as  a boundary  condition  for  this  Sn 
calculation.  This  S„  calculation  process  is 
continued  until  the  standard  Sn  convergence 
criteria  are  meet.  The  converged  Sn  flux 
solution  from  1=1  to  I=IMMP  (®Sr'n  i-ri)  is 


The  final  flux  generated  from  this 
hybrid  model  calculation  is  composed  of  two 
kinds  of  solution  fluxes,  as  in  Approach  1. 
These  include  the  Sn  approximation  flux  from 
1=1  to  I=IMMP  (0  s",  1|M>))  generated  in  step  b, 
and  the  diffusion  theory  flux  from  I=IMMP  to 
I=IM  (dD(|W  ,HJ) , generated  in  step  a,  as  shown 

In  another  possible  approach  to  reach  the 
goal  of  hybrid  model  implementation,  the  Sn  and 
diffusion  approximation  calculations  are 
performed  simultaneously  in  each  iteration  as 
shown  in  Figure  4-3.  In  order  to  keep  the 
continuity  of  flux  at  the  physical  and 


mathematical  boundaries,  it  is  necessary  to  define 
suitable  boundary  conditions.  At  the  physical 
interface,  the  normal  boundary  conditions  for  the 
continuity  of  angular  flux  is  to  be  used,  i.e., 

(ze,E,d)  - <t>2  (f,,B,Q) 


where  r,  is  the  position  of  the  physical  interface 
(see  Figure  4-3) . 

At  the  mathematical  interface  the  neutron 
flux  will  be  essentially  linearly  anisotropic  to 
isotropic  (see  Chapter  3 for  more  details) . 
Therefore,  at  this  interface  another  set  of 
boundary  conditions  could  be  used  to  maintain  the 
continuity  such  as 


J«„  <?„)  - J Sit  (r„> 


j;„  (r„>  - Jmi  lzm)  <4'3) 


l>(f„,Q)-es-^  WiUjLt.fJ  (4.4) 

l*i<o 

Ji«y  -f^daa^(rm,a)  e.-T  {4_5) 

l*J>0 

Jdl/(x'„)  - --i  <Mr„)  T («-«> 

e.-J(iO  - (J*(.f)  - 1r(i’)  «-7) 

- vft  flw(r,  Q)  dO  (4-8) 

(4-9) 


is  the  position  of  the  mathematical  interface, 

S„  approximation  for  the  partial  current  in  the 
negative  direction, 

Sn  approximation  for  the  partial  current  in  the 
positive  direction, 

diffusion  theory  approximation  for  the  partial 
current  in  the  negative  direction, 

diffusion  theory  approximation  for  the 
current  in  the  positive  direction. 


interface  ( i . e . , 


interface  between 


and  SjP,  approximations) , the  normal  boundary 
condition  for  the  continuity  of  angular  flux 
[Equation  (4-1)]  is  to  be  used.  At  the 
mathematical  interface  (the  interface  between  the 
S2P,  and  diffusion  approximations)  another  set  of 
boundary  conditions,  such  as  Equations  (4-2)  and 
(4-3),  should  be  used  to  maintain  the  continuity. 

The  disadvantage  of  both  Approachs  3 and  4 is 
that  the  diffusion  theory  calculation  is  forced  to 
undergo  the  same  number  of  outer  iterations  as  the 
S„  approximation;  as  indicated  in  Chapter  3, 
diffusion  theory  generally  requires  few  iterations 
to  converge  than  does  an  Sn  calculation. 

One  pjmensjonal  Sn-AlbedO-«BflS]. 

Another  possible  approach  to  accelerate  one-dimensional 
S„  transport  theory  calculations  in  the  GCR  system  is  to  use 
a one-region,  one-dimensional  Sn  transport  theory 
approximation  to  calculate  the  flux  within  the  gas  core;  an 
albedo  boundary  condition  is  then  used  to  treat  the  effects 
of  the  reflector  on  the  core. 

The  advantage  of  this  approach  is  that  the  computation 
time  of  the  S„  transport  theory  calculation  should  be  short 
since  calculations  are  performed  only  in  the  gas  core 


This  approach  can  be  used  to  perform  k,„  or  critical 
mass  calculations  and  to  obtain  estimates  of  the  flux 
distribution  within  the  gas  core  region.  However,  this 
method  gives  no  information  directly  about  the  flux  behavior 
in  the  reflector,  although  such  information  may  be  obtained, 
if  desired,  by  solving  additional  equations  at  the  cost  of 
additional  computer  time. 

This  approach  appears  to  be  promising,  but  in  general 
it  will  likely  have  a somewhat  larger  error  associated  with 
it  than  the  "hybrid"  method  explained  in  Approach  1; 
however,  it  should  also  be  faster  than  the  "hybrid" 
approach . 

Generalized  neutron  Diffusion  Equation 

Many  techniques  have  been  suggested  to  accelerate  S„ 
transport  approximation  calculations.  One  approach  is  to 
obtain  a generalized  neutron  diffusion  equation,  the 
solution  of  which  is  an  accurate  approximation  to  the 
transport  solution.  A number  of  methods  have  as  their  aim 
the  improvement  of  the  diffusion  coefficient  of  the  neutron 
diffusion  equation.  This  improvement  of  the  diffusion 
coefficient  leads  to  closer  agreement  between  the  diffusion 
theory  calculated  results  and  the  actual  transport  results. 
Therefore,  assuming  that  the  computation  time  for  solving 
the  new  diffusion  equation  remains  relatively  unchanged  and 
that  the  solution  is  transport  accurate,  the  computational 
time  savings  can  be  substantial  over  that  necessary  for  the 


full  transport  solution. 


improved  diffusion 


coefficient  is  evaluated  from  transport  information 
superior  to  that  used  in  conventional  diffusion  theory. 

Some  of  these  methods  are  outlined  below: 

1.  One  technique  is  to  generate  a diffusion 
coefficient  which  depends  upon  direction  but  in 
the  form  of  a tensor  which  is  evaluated  from 
integral  transport  kernels  [45].  However,  the 
diffusion  tensor  is  full,  resulting  in  a more 
complicated  coupling  of  derivatives  of  the  flux 
than  in  the  conventional  diffusion  equation  and 
the  transport  kernels  are  easily  utilized  only  in 
geometrically  simple,  homogeneous  media. 
Therefore,  this  approach  presents  no  promise  of 
being  able  to  provide  a one  to  two  order  of 
magnitude  improvement  in  computation  time  for  GCR 
transport  theory  calculations. 

2.  A second  technique  is  to  evaluate  the 
diffusion  coefficient  from  auxiliary  transport 
information.  This  information  is  determined  from 
a full  transport  solution  for  a reference 
configuration  of  a system  for  which  a series  of 
calculations  of  perturbed  configurations  is  to  be 
performed.  In  one  case  this  information  is 
utilized  in  a synthesis  technique  in  which  the 
angular  dependence  is  eliminated;  in  another  case 
the  information  is  used  to  evaluate  a full 


diffusion  tensor  whose  form  has  been  derived  from 
a first  angular  moment  balance  of  the  transport 
equation  [76,77].  This  technique  seems  to  be  very 
complicated  and  time  consuming.  Therefore,  this 
approach,  like  the  above  approach,  is  given  no 
further  consideration. 

A third  technique  is  to  formulate  a 
generalized  diffusion  equation  capable  of  yielding 
an  accurate  solution  for  the  transport  scalar  flux 
[78].  In  contrast  to  the  two  methods  mentioned 
above,  this  equation  is  based  on  a diagonal 
diffusion  tensor  and,  hence,  can  be  solved  by 
conventional  numerical  methods  using  standard 
codes.  The  diagonal  diffusion  tensor  incorporates 
information  from  computationally  cheap  auxiliary 
transport  calculations.  This  method  could  be 
useful  for  two-dimensional  problems  but  the  scope 
of  the  work  being  performed  as  part  of  this 
research  does  not  extend  beyond  one  dimension. 

In  this  formulation  the  diffusion  equation  is 
generalized  to 

-V-B(r,E)-V®(r,B)*S,(r,B)®(r,B)-S0(r,B)  (4-10) 


a diagonal  diffusion  t 


components 


-i-  fa  ,¥  (r,Q,  E)  dQ 
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*(r,n,E)  = auxiliary  transport  flux 


t0(r,E)  = <r(r,n,E)  dn  is  the  transport  scalar 


£R(r,E)  = removal  cross  section 
e(r,E)  - diffusion  scalar  flux 

In  the  above  expressions,  i and  j represent 
the  directions  of  the  channels.  A channel  is  that 
portion  of  the  two-dimensional  system  enclosed  by 
lines  parallel  to  one  of  the  coordinate 
directions . 

The  numerical  form  of  this  generalization  is 
sufficient  to  make  the  diffusion  equation 
compatible  with  the  transport  equation,  in  the 
sense  that  the  transport  scalar  flux  can  be 
obtained  as  a solution  to  the  generalized 
diffusion  equation. 

The  numerical,  generalized  diffusion  equation 
is  derived  from  the  numerical  Sn  equation  by 
defining  an  approximating  function  for  the  angular 
flux  which  is  structured  to  yield  the  diffusion 
equation  form.  This  approximating  function  is 
substituted  into  the  Sn  equations,  and  summed  over 


all  angles  to  obtain  the  desired  balance  equation. 
This  derivation  of  the  diffusion  equation  is 
superior  to  the  alternative  of  postulating  a 
generalized  Fick's  law  because  we  can  show  the 
relationship  between  the  scalar  flux  solution  of 
the  transport  equation  and  the  flux  solution  of 
the  diffusion  equation.  The  conditions  that  make 
the  generalized  diffusion  solution  an  accurate 
approximation  of  the  transport  scalar  flux  can 
then  be  determined.  The  auxiliary  transport  flux 
can  be  generated  and  improved.  The  following 
steps  describe  the  procedure  used  in  two- 
dimensional  calculations  to  compute  the  auxiliary 
transport  fluxes  and  give  a plausible  argument  to 
indicate  that  the  auxiliary  information  is  indeed 
improved: 

a.  Compute  the  conventional,  two-dimensional 
diffusion  flux. 

b.  Compute  one-dimensional  transport  fluxes  for 
each  channel  containing  a dominant  source, 
using  the  previous  conventional  flux  to 
estimate  the  transverse  leakage. 

c.  Compute  the  two-dimensional,  generalized 
diffusion  flux  with  the  diffusion 
coefficients  evaluated  from  the  above 
auxiliary  transport  information  in  those 
regions  and  directions  of  the  system 


containing  a channel,  and  from  the 
conventional  diffusion  coefficient  for  the 
remaining  regions  and  directions. 

d.  Compute  a two-dimensional  estimate  of  the 
transport  flux  using  as  a given  source  the 
isotropic  scattering  and  fission  rates 
calculated  from  the  previous  generalized 
diffusion  solution. 

e.  Compute  another  two-dimensional  generalized 
diffusion  flux  with  the  diffusion  coefficient 
evaluated  from  the  two-dimensional  auxiliary 
transport  information  in  step  d. 

Steps  d and  e may  repeated,  if  necessary,  for 
accuracy. 

In  step  b the  one-dimensional  transport  flux 
calculation  is  restricted  to  channels  of  the  system 
containing  a source  region.  The  flux  computed  for  such  a 
channel  is,  therefore,  a representation  of  the  transport 
flux  for  that  portion  of  the  two-dimensional  system  and  that 
coordinate  direction  included  within  the  channel . 

It  is  possible  to  cover  the  entire  two-dimensional 
system  with  one-dimensional  flux  estimates  by  choosing  a 
series  of  channels  to  encompass  the  system  in  each  of  the 
coordinate  directions. 

The  first  step  of  the  above  five  steps  is  to  compute 
the  conventional,  two-dimensional  diffusion  flux. 

Therefore,  this  approach  is  inappropriate  for  the  GCR 


system,  because,  in  general,  diffusion  theory  is  invalid  in 
the  gas  core  region  (see  chapter  3) . In  order  to  make  this 
approach  appropriate  for  two-dimensional  GCR  system 
calculations,  the  following  modified  set  of  steps  might  be 
employed: 


1.  compute  an  estimate  of  the  transport  flux  in 
the  gas  core  region  in  the  primary  direction  from 
a simple  one-dimensional  Sn  calculation  in  which 
a very  low  level  of  convergence  is  imposed  ( e . g . , 

2.  Use  an  albedo  approach  to  estimate  the  core 
flux  in  the  other  direction;  1-D  Sn  transport 
theory  is  used  in  the  core. 

3.  Use  the  fluxes  from  Steps  1 and  2 to  evaluate 
the  core-average  diffusion  coefficient.  This 
diffusion  coefficient  cannot  be  determined  unless 
there  is  a reasonable  degree  of  grayness  yielding 
a definite  curvature  of  the  flux  inside  the  gas 
core.  If  there  is  no  such  flux  curvature  inside 
the  gas  core,  a local  diffusion  coefficient  could 
be  determined  for  each  point  inside  the  gas  core. 
However,  determination  of  a local  diffusion 
coefficient  will  require  a relatively  long 
computation  time.  Hence,  from  a practical 
viewpoint  the  condition  needed  to  employ  this 


approach  is  to  have  a reasonable  degree  of 
grayness  in  the  core. 

4.  use  the  generalized  diffusion  eguation  (and 
the  diffusion  coefficient  evaluated  in  step  3)  to 
compute  the  diffusion  scalar  flux  in  two 
dimensions. 

The  diffusion  flux  calculated  in  Step  4 is  an 
improved  diffusion  scalar  flux  because  the 
modified  diffusion  coefficient  (Step  3)  is 
obtained  from  Sn  transport  theory  fluxes. 

5.  After  all  iterations,  the  generalized 
diffusion  equations  should  yield  a transport 
accurate  flux.  This  is  because  the  improved 
diffusion  coefficient  is  evaluated  from 
transport  information. 

The  SK,  Method 

The  SK,  method  is  another  possible  approach  which  might 
be  used  to  shorten  the  calculation  time  for  GCR  transport 
theory  calculations  [79].  The  SK,  method  is  a high-order 
transport  approximation  introduced  for  solving  the  integral 
transport  eguation  [80-82] . The  method  relies  on 
approximating  the  integral  transport  kernels  by  a sum  of 
diffusion-like  kernels.  The  integral  equation  is  equivalent 
to  a set  of  coupled  second-order  differential  equations  for 
which  blackbody  and  reflecting  boundaries  are  established. 
The  original  formulation  of  the  SK,  approximation  has  a 


defect  when  applied  to  heterogeneous  problems.  Therefore, 
the  SK,  approximation  is  slightly  altered  to  solve  the 
integral  transport  equation  for  heterogeneous  systems  [83]. 
This  technique  can  also  be  applied  to  problems  with  P, 
scattering.  As  introduced  in  references  [79,83],  the  SK, 
method  was  tested  against  a variety  of  benchmark  problems. 
It  is  noticed,  from  these  applications,  that  some  of  the 
benchmark  problems  for  which  this  method  proved  effective 
have  properties  like  those  of  gas  core  reactors.  For 
example,  in  one  of  the  benchmark  problems  the  ratio  (Ej/S,) 
is  0.25.  This  ratio  is  comparable  to  that  of  typical  GCR's 
(see  Table  3-1) . This  method  may,  therefore,  be  applicable 
to  the  GCR.  The  calculations  in  references  [72,76]  were 
performed  on  an  IBM  PC-XT-compatible  machine  that  was 
equipped  with  a mathematical  coprocessor.  The  execution 
time,  for  one  of  these  benchmark  problems,  with  spherical 
geometry,  for  SKy  orders  2,  3,  4,  and  5 are  27  seconds,  46 
seconds,  73  seconds,  and  115  seconds,  respectively.  The 
execution  times  for  S12  calculations  with  ANISIN,  for  the 
same  benchmark  problem,  are  4.04  minutes,  7.35  minutes,  and 
4.20  minutes  for  slab,  cylindrical,  and  spherical 
geometries,  respectively.  Comparison  of  these  execution 
times  shows  that  the  SK,  method  can  be  much  faster  than  the 
standard  sn  calculation.  Also,  it  is  shown  for  these 
benchmark  problems  that  the  SKy  method  results  are  more 
than  the  Sn  results.  Therefore,  it  is  possible 


that  the  SK,  method  might  be  superior  to  the  S„ 
approximation  for  GCR  calculations. 

SK,, -Diffusion  Hybrid  Model 

This  hybrid  model  depends  on  the  use  of  an  initial 
diffusion  theory  calculation  for  both  the  gas  core  and 
reflector.  Then  the  SK,  calculation  is  used  for  the  gas 
core  region  and,  if  necessary,  for  the  first  few  mean  free 
path  inside  the  reflector.  The  steps  required  for  this  SK,- 
diffusion  hybrid  model  are  similar  to  the  steps  explained  in 
the  first  approach  of  the  Sn-diffusion  hybrid  model  (i.e., 
similar  to  the  approach  used  in  this  dissertation) . This 
hybrid  model  could  be  faster  than  the  S„-diffusion  hybrid 
model  because  the  SK,  method  is  much  faster  than  the  Sn 
approximation  calculation  as  explained  in  the  above  section. 

Sce^Biiaensicnal  SH,rAibedo  Model 

Another  possible  approach  to  accelerate  the  standard 
one-dimensional  S„  transport  theory  calculation  in  the  GCR 
system  is  to  use  a one-dimensional  SK,  calculation  within 
the  gas  core  with  an  albedo  boundary  condition  used  to  treat 
the  effect  of  the  reflector  on  the  core. 

The  advantage  of  this  approach  is  that  the  computation 
time  of  the  SK,  calculation  should  be  very  short  since 
calculations  are  performed  only  in  the  gas  core  region. 


approach  appears 


promising. 


a somewhat  larger  error  associated  with  it  than  the  SK,- 
dif fusion  hybrid  model;  however,  it  should  also  be  faster 
than  the  hybrid  model. 

On  the  other  hand,  this  approach  can  be  used 
successfully  for  many  different  purposes.  For  example,  it 
can  be  used  to  perform  K,()  or  critical  mass  calculations  and 
perform  flux  calculations  within  the  gas  core  region. 
However,  this  method  gives  no  information  directly  about  the 
flux  behavior  in  the  reflector. 

Vectorization  and  MultitasKing 

Vectorization  and/or  Multitasking  are  two  techniques 
which  may  be  used  to  accelerate  FORTRAN  codes  or  to  improve 
the  acceleration  of  any  third  technique. 

Vector  processing  is  a technique  for  reducing  the 
amount  of  processor  time  required  to  run  a program.  It  can 
be  applied  to  a FORTRAN  code  segment  such  as  a DO-LOOP  that 
performs  the  same  sequence  of  operations  on  successive 
elements  of  arrays.  Vector  processing  reduces  the  processor 
time  required  to  execute  such  loops  by  overlapping,  or 
pipelining,  the  processing  for  the  array  elements. 

The  VS  FORTRAN  Version  2 compiler  can  automatically 
generate  vector  instructions  in  appropriate  places  in  the 
object  module  as  it  compiles  the  VS  FORTRAN  (language  level 
66  or  77)  source  program;  this  process  is  called 


Vectorization.  However,  this  compiler  examines  each  DO-LOOP 
and  selects  for  vectorization  only  those  statements  that  can 
be  safely  and  effectively  transformed.  Statements  that 
cannot  be  safely  vectorized,  or  would  not  run  faster  in 
vector  mode,  are  compiled  into  standard  scalar  instructions. 

Multitasking  is  designed  to  improve  turnaround  time  for 
computer-intensive  applications.  The  amount  of  turnaround 
improvement  that  can  be  achieved  depends  on  the  application. 
The  primary  factors  affecting  turnaround  performance  are  the 
number  of  processors  and  the  percentage  of  the  application 
processing  that  can  be  performed  in  parallel.  It  can  be 
used  on  programs  that  can  be  divided  into  two  or  more 
independent  subtasks  and  run  on  a multiprocessor 
configuration.  This  allows  the  subtasks  to  be  executed 
concurrently,  when  the  system  workload  permits,  thereby 
potentially  reducing  the  turnaround  time  required  by  the  job 
(see  Appendix  c for  more  information  about  Vectorization  and 


Multitasking  techniques) . 
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theory  flux  calculations  required  when  a fission  density 
guess  is  supplied,  or  when  flagged  in  the  input.  The 
standard  DT  subroutine  generates  scalar  flux  values  at  the 
center  of  each  interval  and  cannot  calculate  partial 
currents.  Therefore,  this  subroutine  was  modified  to 
calculate  the  value  of  the  scalar  flux  at  the  right  boundary 
of  the  mathematical  interface  interval  (i.e.,  at  I=IMMP)  and 
the  partial  currents  (J-  and  J‘)  for  each  energy  group  at 
the  right  boundary  of  this  interval.  These  partial  currents 
values  are  then  used  in  the  OUTERS  subroutine  to  calculate 
the  group  dependent  albedos  which  are  used  as  the  outer 
boundary  conditions  for  the  Sn  portion  of  the  hybrid 
calculation.  The  scalar  fluxes  generated  in  the  diffusion 
portion  of  the  hybrid  calculation  range  from  1=1  to  I=IM. 
However,  initial  flux  guess  values  for  the  Sn  portion  of  the 
hybrid  calculation  are  only  needed  from  1=1  to  I=IMMP.  The 
subroutine  WRTF  which  write  out  "restart"  fluxes  was 
therefore  modified  so  that  the  converged  diffusion  theory 
fluxes  which  are  used  as  restart  or  initial  guess  flux 
values  for  the  Sn  calculation  range  only  from  1=1  to  I=IMMP. 
The  discrete  ordinates  angular  flux  calculations  are 
performed  in  the  INNER  subroutine.  It  uses  the  current 
scalar  flux  and  flux  moments  to  calculate  the  scattering 
source  for  the  current  group.  INNER  then  loops  over  the 
angels,  calculating  the  angular  flux  at  each  mesh  boundary 
using  a weighted  diamond  difference  scheme,  traversing  the 
mesh  in  the  direction  of  the  current  angle.  Figure  5-1 


shows  the  general  flow  diagram  for  hybrid  calculations 
performed  using  the  modified  XSDRNPM-S  computer  code. 

In  order  to  make  the  structure  and  the  implementation 
of  the  hybrid  model  clear  and  understandable,  the  flow 
diagrams  of  the  main  three  subroutines  involved  in  this 
hybrid  calculation  are  shown  in  Appendices  D,  E,  and  P. 
Appendix  D includes  the  flow  diagram  of  the  OUTERS 
subroutine,  while  Appendix  E shows  the  flow  diagram  of  the 
DT  subroutine,  and  Appendix  F includes  the  flow  diagram  of 
the  INNER  subroutine. 

Neutronic  Analysis  of  a Solid  Core  Reactor 


Neutronics  calculations  have  been  performed  on  a solid 
core  reactor  with  an  H20  moderator/reflector  (see  Figure 
1-2) . This  benchmark  problem  was  selected  for  this 
particular  study  because  both  S„  transport  theory  and 
diffusion  theory  approximations  can  be  successfully  applied 
to  this  problem,  and  give  essentially  the  same  results. 
Neutronics  calculations  performed  on  the  solid  core  reactor 
in  one  dimension  using  these  approximations  have  verified 
that  the  error  associated  with  the  diffusion  theory 
approximation  as  compared  to  Sn  transport  theory 
calculations  is  small  for  this  benchmark  problem.  Hence, 
the  "hybrid"  model  (the  model  developed  in  this  work  to 
obtain  reduced  computation  times  relative  to  the  pure  S 
approximation  for  GCR  calculations)  should  also  give  the 


same  results  for  this  system.  This  solid  core  benchmark 
problem  is  to  serve  as  a diagnostic  tool  to  verify  whether 
the  "hybrid"  neutronics  model  that  will  be  used  in  the 
system  of  interest,  i.e.  the  gas  core  reactor,  has  been 
properly  implemented. 

In  this  study,  NITAWL  is  used  to  generate  a 123  fine 
group  neutron  cross  section  library  for  this  particular 
problem.  The  123  fine  group  neutron  cross-sections  (30 
thermal  groups)  are  collapsed  down  to  a 26  energy-group 
structure  (6  thermal  groups)  and  then  to  a 4 energy-group 
structure  (one  thermal  group)  with  XSDRNPM  [2]. 

XSDRNPM  has  the  ability  to  perform  both  Sn  transport 
theory  and  diffusion  theory  calculations.  It  is  used  to 
perform  both  one-dimensional  s2P,,  S,P5,  SsP3,  and  diffusion 
theory  calculations  on  the  solid  core  reactor  benchmark 
problem.  For  the  benchmark  problem  with  4 energy  group 
neutron  cross  sections,  the  diffusion  theory  approximation 
is  found  to  overpredict  k,„  by  less  than  1%  relative  to  the 
SjP,,  s,P3,  and  S6P3  results  in  the  slab  geometry.  In 
addition,  26  energy  group  calculations  have  been  performed 
using  both  the  diffusion  theory  and  an  StP5  approximation 
for  slab  geometry.  The  diffusion  theory  approximation,  with 
26  energy  groups,  is  also  found  to  overpredict  kr((  by  less 
than  1%  relative  to  the  S6P3  results.  These  calculations 
have  been  repeated  for  different  dimensions  and  geometries. 
Tables  5-1,  5-2,  and  5-3,  show  the  results  of  these 
calculations  for  three  different  H;Q  reflector  thicknesses. 
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The  comparison  of  neutron  multiplication  factor  (k^,) 
values,  as  shown  in  these  tables,  indicates  that  for  this 
benchmark  problem  diffusion  theory  provides  good  accuracy 
(<1%  error  in  ke()  for  slab  geometry) . An  examination  of  the 
total  number  of  iterations  and  the  CPU  computation  times  in 
these  tables  shows  that  the  diffusion  theory  approximation 
is  much  faster  than  the  S„  transport  theory  calculation  (by 
a factor  of  up  to  10  times  for  4 groups  and  a factor  of  up 
to  about  8 times  for  26  groups) . Such  savings  in 
computation  times  are  much  more  significant  in  the  case  of 
multigroup  calculations.  For  example,  from  Table  5-3  the 
diffusion  theory  calculation  is  10  times  faster  for  four 
groups  than  the  S6PS  calculation  and  the  saving  in  CPU  time 
is  about  7 seconds.  On  the  other  hand,  for  the  26  group 
problem  although  the  diffusion  theory  calculation  is  only 
about  8 times  faster  than  the  S6P5  calculations,  the  saving 
in  CPU  time  is  almost  16  seconds. 

XSDRNPM-S  has  also  been  used  to  perform  S„  and 
diffusion  theory  calculations  on  this  benchmark  problem. 
These  calculations  were  performed  with  4 energy  group 
neutron  cross  sections.  The  diffusion  theory  approximation 
is  again  found  to  overpredict  kffff  by  less  than  1%  relative 
to  the  SjP,,  S4Pj,  and  S6P3  results  in  slab  geometry. 

XSDRNPM-S  has  also  been  used  to  perform  s6P5  and  diffusion 
theory  calculations  using  26  energy  groups.  The  diffusion 
theory  approximation,  with  26  energy  groups,  is  also  found 
to  overpredict  ke(,  by  less  than  1%  relative  to  the  s6P3 


results.  Table  5-4  shows  the  results  of  the  Solid  Core 
Reactor  benchmark  calculations  performed  with  the  XSDRNPM-S 
computer  code.  As  expected,  the  results  of  the  calculations 
performed  by  XSDRNPM-S  are  comparable  to  the  results  of  its 
predecessor,  XSDRNPM. 

Neutronic  calculations  are  also  performed  on  benchmark 
problem  using  the  hybrid  model  developed  in  this  work.  The 
kI((  results  from  the  hybrid  model  and  the  standard  S„  and 
diffusion  theory  approximations  are  compared  with  one 
another  in  Table  5-5.  This  table  shows  that  the  hybrid 
model  is  found  to  underpredict  ke),  only  about  0.01%  relative 
to  the  standard  Sn  approximation.  An  examination  of  the 
total  number  of  iterations  and  the  CPU  calculation  times  in 
Table  5-5  shows  that  the  hybrid  model  is  roughly  2.6  times 
faster  than  the  standard  Sn  approximation  calculation  for 
this  benchmark  problem  configuration.  These  results 
indicate  that  the  hybrid  model  developed  in  this  work  has 
been  properly  implemented  in  XSDRNPM-S. 
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CHAPTER 


R NEUTRONICS  ! 


UDIES  USING  THE  XSDRNPM-S  COMPUTER  C 
E DEVELOPED  HYBRID  MODEL 


Iteutiwils  Studies  of  the  SCR  Vail 


The  XSDRNPM-S  computer  code  solves  either  the  S„ 
approximation  to  the  transport  equation  or  the  diffusion 


performed  on  the  GCR  using  the  XSDRNPM-S  computer  code. 


configurations,  the  diffusion  theory  approximation  is  found 
to  overpredict  ke((  by  less  than  1%  relative  to  an  StPs 
approximation.  An  examination  of  the  total  number  of 


groups) 


total  savings  in  computation  time  are  much  more  significant 
in  the  case  of  multigroup  calculations  and/or  for  thick 
reflectors.  The  diffusion  theory  calculations  for  the  gas 
core  reactor  were  found  to  converge  in  just  a few  iterations 
(only  two  or  at  most  three  iterations  are  generally  needed 
to  reach  essentially  total  convergence) . In  contrast,  the 
Sn  approximation  calculations  generally  required  a 
significant  numbers  of  total  iterations  (i.e.,  greater  than 
40)  to  reach  a convergence  level  of  about  10'4.  For  both 
the  diffusion  theory  and  S„  calculations,  a flat  fission 
density  guess  was  used  to  start  the  calculations.  Because 
the  diffusion  theory  calculation  converges  very  rapidly  in 
this  case,  comparison  of  the  diffusion  theory  computation 
time  with  that  of  the  corresponding  Sn  calculation  can  be 
somewhat  misleading.  Therefore,  a bad  initial  flux  guess  is 
used  for  the  diffusion  calculations  in  order  to  increase  the 
number  of  iterations  and  calculation  time.  At  the  same 
time,  in  order  to  fairly  compare  the  Sn  results  with  the 
diffusion  theory  results,  the  same  bad  guess  is  then  also 
used  for  the  corresponding  Sn  calculation. 

The  results  of  the  comparative  calculations  performed 
with  the  XSDKNPM-S  computer  code  are  comparable  to  the 
results  generated  with  its  predecessor,  the  XSDRNPM  computer 


comparison) . 


core,  externally  moderated  “cavity"  reactor  systems  fueled 
with  OF,  (see  Figure  1-1)  using  the  hybrid  model  developed 
in  this  work.  The  general  description  of  the  hybrid  model 
and  of  how  it  performs  hybrid  calculations  are  presented  in 


Chapters  4 and  5.  Different  GCR  problems  are  examined  using 


results  from  the  hybrid  model  are  compared  with  those  of  the 
standard  S„  and  diffusion  theory  approximations. 


tested  for  many  different  mathematical  interface  positions 


cm,  6 cm,  4 cm,  and  2 cm  inside  the  reflector) . The  ke(( 


diffusion  theory  approximations  are  compared  with  one 
another  in  Table  6-2.  This  table  shows  that  the  hybrid 
model  is  found  to  overpredict  ke(,  by  less  than  1%  relativ 
to  the  standard  S„  approximation  when  the  mathematical 


cm  inside  the 


hybrid  model  underpredicts  ke((  by  less 


sflectc 


The  flux  has  been  normalized  to  1 (fission/sl  and 
multiplied  by  104. 


Figure  6-1  The  flux  behavior  from  the  hybrid  model  com- 
pared with  the  fluxes  from  both  Sn  and 
diffusion  theory  approximations  for  the  case 
when  the  mathematical  interface  is  located 
14  cm  inside  the  reflector  and  D.F.=5.0. 


- The  flux  has  been  normalized  to  1 (fission/s)  and 

multiplied  by  104. 


Figure  6-2  The  flux  behavior  from  the  hybrid  model  com- 
pared with  the  flux  from  the  Sn  approximation 
for  the  case  when  the  mathematical  interface 
is  located  14  cm  inside  the  reflector  and 


up  to  I=IM) . It  is  also  found  that  these  two  flux  solutions 
exhibit  even  better  continuity  when  the  hybrid  flux  is  used 
from  1=1  up  to  I=IMMP-1  and  the  diffusion  theory  flux  is 
used  from  I=IMMP-1  up  to  I=IM  as  shown  in  Figure  6-3.  When 
IMMP  is  located  at  12  cm  and  10  cm  inside  the  reflector,  the 
hybrid  model  gives  the  same  good  flux  behavior  as  for  the  14 
cm  case  discussed  above.  Figures  6-4  to  6-9  show  these 

For  the  case  when  IMMP  is  located  at  8 cm  inside  the 
reflector,  the  results  of  the  hybrid  model  are  still  good 
but  start  to  deviate  slightly  from  that  of  the  standard  S„ 
approximations  as  shown  in  Figures  6-10  to  6-13.  Figure  6- 
10  shows  the  flux  solutions  from  the  hybrid  model,  Sn,  and 
diffusion  theory  approximations  together.  The  flux 
solutions  from  only  the  hybrid  model  and  S„  approximation 
are  shown  in  Figure  6-11  while  the  flux  solution  from  only 
the  hybrid  model  and  diffusion  theory  are  shown  in  Figure 
6-12.  For  this  case  also,  the  continuity  between  the  hybrid 
model  flux  and  that  from  diffusion  theory  is  better  when  the 
transition  is  taken  at  I=IMMP-1  ( i . e . , when  the  hybrid  model 
flux  solution  is  used  from  1=1  up  to  I»IMMP-1  and  the 
diffusion  theory  flux  solution  is  used  from  I=IMMP-1  up  to 
I=IM;  see  Figure  6-13) . 

mathematical  interface  is  located  at  6 cm  or  less  inside  the 
reflector,  for  this  particular  configuration,  is  related  to 
the  nonlinear  anisotropic  behavior  of  the  flux.  As 


The  flux  has  been  normalized  to  1 (fission/s)  and 

multiplied  by  104. 


Figure  6-3  The  flux  from  the  hybrid  model  from  1=1  up 
to  1 = (IMMP-1)  merged  together  with  the 
diffusion  theory  flux  from  I B (IMMP-1)  up 
to  I = IM  for  the  case  when  the  mathematical 
interface  is  located  14  cm  inside  the 
reflector  and  D.F.=5.0. 
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The  £lux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10^. 


The  flux  behavior  from  the  hybrid  model  com- 
pared with  fluxes  from  both  Ss  and  diffusion 
theory  approximations  for  the  case  when  the 
mathematical  interface  is  located  12  cm  inside 
the  reflector  and  D.F.  - 5.0. 


- The  flax  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  104. 


Figure  6-5  The  flux  behavior  from  the  hybrid  model  com- 
pared with  the  flux  from  the  Sn  approximation 
for  the  case  when  the  mathematical  interface 
is  located  12  cm  inside  the  reflector  and 
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The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10*. 


Figure  6-6  The  flux  from  the  hybrid  model  from  I = 1 up 
to  I = (IMMP-1)  merged  together  with  the 
diffusion  theory  flux  from  I = (IMMP-1)  up  to 
I = IM  for  the  case  when  the  mathematical 
interface  is  located  12  cm  inside  the  re- 


- The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10*. 

Figure  6-7  The  flux  behavior  from  the  hybrid  model  com- 
pared with  fluxes  from  both  Sn  and  diffusion 
theory  approximations  for  the  case  when  the 
mathematical  interface  is  located  10  cm  in- 
side the  reflector  and  D.F.  = 5.0. 
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- The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10^. 


Figure  6-8  The  flux  behavior  from  the  hybrid  model  com- 
pared with  the  flux  from  the  Sn  approximation 
for  the  case  when  the  mathematical  interface 
is  located  10  cm  inside  the  reflector  and 


- The  flux  has  been  normalised  to  1 (fission/s)  and 
multiplied  by  104. 


Figure  6-9  The  flux  from  the  hybrid  model  from  1*1  up 
to  I = tIMMP-1)  merged  together  with  the 
diffusion  theory  flux  from  1 = (IMMP-1)  up 
to  1 = IM  for  the  case  when  the  mathematical 
interface  is  located  10  cm  inside  the  re- 


2.0 


Radial  Position  (cm) 


The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10^. 


Figure  6-10  The  flux  behavior  from  the  hybrid  model  com- 
pared with  the  fluxes  from  both  Sn  and 
diffusion  theory  approximation  for  the  case 
when  the  mathematical  interface  is  located 
8 cm  inside  the  reflector  and  D.F.-5.0. 


- The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  104 . 

Figure  6-11  The  flux  behavior  from  the  hybrid  model  com- 
pared with  the  flux  from  the  Sn  approxima- 
tion for  the  case  when  the  mathematical 
interface  is  located  8 cm  inside  the  reflec- 


The  flux  has  been  normalized  to  1 (fission/s)  and 

multiplied  by  104. 


Figure  6-12  The  flux  behavior  from  the  hybrid  model  com- 
pared with  the  flux  from  diffusion  theory 
for  the  case  when  the  mathematical  interface 
is  located  8 cm  inside  the  reflector  and 


- The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10^. 


Figure  6-13  The  flux  from  the  hybrid  model  from  1 = 1 
up  to  1 = (IMMP-1)  merged  together  with 
the  diffusion  theory  flux  from  I = IMMP-1 
up  to  I = IM  when  th4e  mathematical  inter- 
face is  located  B cm  inside  the  reflector 


discussed  in  chapter  3,  a certain  degree  of  nonlinear 
anisotropic  behavior  can  be  exhibited  by  the  core  flux  in  a 
GCR  and  this  nonlinear  anisotropic  behavior  extends  a 
certain  distance  inside  the  reflector.  This  distance 
increase  with  increasing  density  factor  for  the  gas  core. 
This  nonlinear  anisotropic  behavior  of  the  flux  for  the 
first  few  Bean  free  path  inside  the  reflector  renders 
diffusion  theory  invalid  for  the  same  distance  inside  the 
reflector. 

The  albedo  values  used  in  the  hybrid  model  are 
generated  from  the  converged  diffusion  theory  scalar  flux. 

If  diffusion  theory  is  not  valid  at  a certain  location,  it 
will  give  wrong  albedo  values  and  these  wrong  albedo  values 
will  cause  the  hybrid  model  to  give  poor  results.  This  is 
precisely  what  happens  when  IMMP  is  located  at  6 cm  or  less 
inside  the  reflector  for  this  particular  configuration. 

In  order  to  make  the  above  discussion  clearer,  another 
gas  core  configuration  with  D.F.=3.0  instead  of  D.F.-5. 0 is 
examined.  This  density  factor  (i.e.,  D.F.=3.0)  is  used 
because  it  yields  a system  which  exhibits  a weakly  nonlinear 
anisotropic  behavior  (see  chapter  3,  Figures  3-5  and  3-6). 
This  means  that  it  is  expected  that  the  hybrid  model  will 
give  better  results  at  IMMP=6  cm  for  D.F.=3.0  than  for 
D.F.=5.0.  Table  6-3  compares  the  results  for  D. F . =3 . 0 with 
the  corresponding  results  for  D.F.=5.0  when  IMMP  is  located 
at  6 cm  and  4 cm  into  the  reflector.  It  is  found  that  the 
hybrid  model  gives  acceptable  results  (error  in  kelf  of  1.4%) 
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when  IMMP  is  located  6 cm  into  the  reflector  with  D.F.“3.0 
and  unacceptable  results  (error  in  ks„  of  almost  45%)  when 
IMMP  is  located  6 cm  into  the  reflector  with  D.F.=5.0.  The 
flux  behavior  prediction  from  the  hybrid  model  for  D.F.=3.0 
is  compared  with  the  flux  behavior  prediction  from  both  the 
S„  and  diffusion  theory  approximations  when  IMMP  is  located 
6 cm  into  the  reflector  in  Figures  6-14  to  6-17.  Figure  6- 
15  shows  that  the  hybrid  model  flux  behavior  is  in  good 
agreement  with  that  of  the  Sn  approximation.  Also,  it  is 
found  that  there  is  relatively  good  continuity  between  the 
flux  from  the  hybrid  model  at  IMMP  and  the  flux  from 
diffusion  theory  at  IMMP  (see  Figures  6-16  and  6-17) . In 
contrast,  Figure  6-18  shows  that  the  hybrid  model  flux 
behavior  is  not  in  good  agreement  with  that  of  the  Sn 
approximation  when  IMMP  is  located  6 cm  into  the  reflector 
for  the  case  when  D.F.-5.0.  Correspondingly,  the  continuity 
is  not  good  between  the  flux  from  the  hybrid  model  at  IMMP 
and  the  flux  from  diffusion  theory  when  D.F.-5.0  (see  Figure 
6-19) . 

When  the  gas  core  reactor  is  operating  at  power,  it 
will  in  general  have  a nonuniform  gas  density  distribution 
in  the  core  due  to  the  nonuniform  temperature  distribution 
inside  the  gas  core  (the  temperature  close  to  the  center  is 
higher  than  the  temperature  in  a region  close  to  the 
reflector) . Therefore,  gas  core  reactor  configurations  with 
nonuniform  gas  density  distributions  in  the  core  are 
examined  to  determine  the  impact  of  nonuniform  density  on 


The  flux  has  been  normalized  to  1 (fission/s)  and 

multiplied  by  104. 


Figure  6-14  The  flux  behavior  from  the  hybrid  model  com- 
pared with  fluxes  from  both  the  Sn  and 
diffusion  theory  approximations  for  the  case 
when  the  mathematical  interface  is  located 
6 cm  inside  the  reflector  and  D.F.=3.0. 
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The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10’. 


Figure  6-15 


The  flux  behavior  from  the  hybrid  model  com- 
pared with  the  flux  from  the  Sn  approxima- 
tion for  the  case  when  the  mathematical 
interface  is  located  6 cm  inside  the  re- 
flector and  D.F.  - 3.0. 


- The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  104. 


Figure  6-16  The  flux  behavior  from  the  hybrid  model  com- 
pared with  the  flux  from  diffusion  theory 
for  the  case  when  the  mathematical  interface 
is  located  6 cm  inside  the  reflector  and 
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The  flux  has  been  normalised  to  1 (fission/s)  and 
multiplied  by  10^. 


Figure  6-17 


The  flux  of  the  hybrid  model  from  I = 1 up 
to  I = IMMP  merged  together  with  the  diffu- 
sion theory  flux  from  I = IMMP  up  to  I = IM 
when  the  mathematical  interface  is  located 
6 cm  inside  the  reflector  and  P.F.  - 3.0. 


The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10^. 


Figure  6-18 


The  flux  behavior  from  the  hybrid  model 
compared  with  the  flux  from  the  Sn  approxi- 
mation for  the  case  when  the  mathematical 
interface  is  located  6 cm  inside  the  re- 
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The  flux  has  been  normalised  to  1 (fission/s)  and 
multiplied  by  104. 


Figure  6-19  The  flux  behavior  from  the  hybrid  model  from 
I - 1 up  to  1 - IMPP  merged  together  with  the 
diffusion  theory  flux  from  I = IMMP  up  to 
I = IM  when  the  mathematical  interface  is 
located  6 cm  inside  the  reflector  and  D.F.= 


- The  flu*  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10^. 


Figure  6-20  Thermal  neutron  flux  variation  versus  radial 
position  for  the  hybrid  model  compared  with 
fluxes  from  both  Sn  and  diffusion  theory 
approximations,  for  D.F.=10.0  in  the  central 
region  of  the  core  and  D.F.=1.0  in  the  outer 


The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  lO^. 


The  flux  of  the  hybrid  model  from  I - 1 up 
to  the  I = IMMP  merged  with  the  diffusion 
theory  flux  from  1 = IMMP  up  to  I - IM,  for 
D. F. =10 . 0 in  the  central  region  of  the  core 
and  D.F.=1.0  in  the  outer  region  of  the 
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- The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  104. 


Figure  6-22  Thermal  neutron  flux  variation  versus  radial 
position  for  the  hybrid  model  compared  with 
fluxes  from  both  Sn  and  diffusion  theory 
approximations,  for  D.F . = 1.0  in  the  central 
region  of  the  core  and  D.F.  = 10.0  in  the 
outer  region  of  the  core. 
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The  flux  has  been  normalized  to  1 (fission/s)  and 
multiplied  by  10^. 


The  flux  of  the  hybrid  model  from  I = 1 up  to 
the  I * IMMP  merged  with  the  diffusion  theory 
flux  from  I = IMMP  up  to  I = IM  for  D.F . = 

1.0  in  the  central  region  of  the  core  and 
D.F.  = 10.0  in  the  outer  region  of  the  core. 


the  flux  from  the  hybrid  model  from  1=1  up  to  I=IMMP  merged 
with  the  diffusion  theory  flux  from  I=IMMP  up  to  I=IM.  This 


The  k,„  results  from  the  hybrid  model  and  the  standard 
sn  approximation  are  compared  with  each  another  for  these 
two  nonuniform  gas  density  cases  in  Table  6-4.  This  table 
shows  that  the  hybrid  model  is  found  to  overpredict  ke(l  by 


different  reflector  thicknesses  are  used  for  this  purpose 


The  hybrid  model  again  overpredicts  ke((  by  less  than  1% 

hybrid  model  is  roughly  five  times  faster  than  the  standard 
Sn  approximation  for  the  case  of  a 100  cm  reflector 
thickness.  For  a 150  cm  reflector  thickness,  the  hybrid 


hybrid  model  are  essentially  identical  to  t 
standard  sn  approximation. 

The  performance  of  the  hybrid  model  is 
for  two  different  energy  group  structures  ( 


o examined 


groups) . 

It  is  found  that  the  hybrid  model  overpredicts  kt„  by 
less  than  1%  for  both  4 energy  groups  and  21  energy  groups 
(see  Table  6-6) . For  4 energy  groups  the  hybrid  model  is 
found  to  be  roughly  five  times  faster  than  the  standard  S„ 
approximation  and  about  five  and  a half  times  faster  than 
the  standard  sn  approximation.  However,  the  total  saving  in 
computation  time  is  much  more  significant  in  the  case  of  21 
energy  groups,  in  spite  of  the  fact  that  the  acceleration 
for  the  two  cases  are  roughly  equal. 


CHAPTER 


FUNDAMENTAL  CONSIDERATIONS  AND  PRELIMINARY 
ANALYSIS  OF  UTVR  SHIELDING  DESIGN 


The  use  of  nuclear  power  in  space,  whether  from 
reactors  or  radioisotopes,  will  be  accompanied  by  the 
presence  of  nuclear  radiations.  Since  these  radiations- 
alpha  particles,  beta  particles,  gamma  rays,  neutrons,  etc.- 
are  potentially  harmful  to  the  astronaut  crew,  delicate 
electronic  components  on  the  spacecraft,  and  scientific  and 
applications  payloads,  space  system  planners  must  include 
radiation  protection  as  a fundamental  design  factor. 

One  of  the  most  fundamental  principles  of  radiation 
protection  is  that  exposure  should  always  be  kept  to  a 
minimum  in  achieving  a particular  goal  or  in  performing  a 
specific  task;  and  exposure  to  an  established  limit  should 
only  be  permitted  when  absolutely  necessary  [84].  The  use 
of  a nuclear  power  source  in  space  is  accompanied  by  the 
need  for  a radiation  shield  to  protect  the  astronaut  crew, 
payload,  and  radiation-sensitive  spacecraft  equipment  from 
the  damaging  effects  of  radiation.  This  is  especially  true 
for  nuclear  reactor  systems  and  the  copious  quantities  of 
neutrons  and  gamma  rays  that  accompany  the  fission  process. 
The  establishment  of  permissible  radiation  levels  to  protect 


personnel  and/or  sensitive  equipment  is  one  of  the  major 
steps  in  the  design  of  a radiation  shield. 

Space  power  systems  utilizing  reactors  must  be 


optimized  with  respect  to  both  size  and  mass  while  also 
meeting  the  performance  specifications  for  neutron  fluence 
and  gamma  dose  to  protect  the  astronaut  crew  (for  manned 
missions)  and  sensitive  spacecraft  equipment,  like  magnets 
and  electronics  parts,  from  damage.  The  reactor  and  shield 
can  represent  a large  fraction  of  the  system  mass,  and  the 
design  of  each  affects  the  other.  Therefore,  it  is 
necessary  to  pay  close  attention  to  optimizing  the  reactor 
and  shield  together  when  optimizing  the  power  system. 
Minimization  of  the  shield  mass,  which  generally  accounts 
for  20%  or  more  of  the  system  mass,  is  very  important. 
Normally,  a shadow  shield  is  located  between  the  reactor 
core  and  the  user  interface  plane.  A choice  of  lithium 
hydride  combined  with  a layer  of  tungsten  has  been  found  to 
be  a good  low  mass  solution  for  space  reactor  shielding 
[85].  In  addition,  the  shield's  internal  heating  must 
produce  tolerable  temperatures.  Finally,  the  shield 
materials  must  remain  physically  and  chemically  stable  for 
the  lifetime  of  the  vehicle.  The  establishment  of 
permissible  radiation  levels  to  protect  personnel  or 
sensitive  equipment  is  one  of  the  major  steps  in  the  design 
of  a radiation  shield.  Figure  7-1  presents  a composite 
description  of  typical  radiation  damage  thresholds  for  man 
and  also  for  certain  materials,  components,  and  systems  that 
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might  be  used  in  nuclear  powered  spacecraft  [86] . when  such 
materials  and  components  are  exposed  to  nuclear  radiations, 
their  properties  are  frequently  changed  in  a variety  of  ways 
that  cannot  be  readily  anticipated  from  nonradiation 
environment  operating  experience. 

Many  factors  have  an  effect  on  the  geometry  and  mass  of 
the  shield,  like: 

a.  the  size  and  nature  of  the  nuclear  power  source; 

b.  the  type  of  radiation; 

c.  the  general  configuration  of  the  spacecraft; 

d.  the  length  of  the  mission;  and  the 

e.  total  permitted  radiation  dose. 

There  are  a variety  of  shield  configurations  that  might 
be  selected  with  the  reactor-powered  space  vehicle.  For 
example,  a shadow  shield  will  only  partially  surround  the 
nuclear  power  source.  The  area  within  the  shadow  of  the 
shield  is  protected,  while  radiation  escapes  to  space  from 
the  reactor's  nonshielded  surfaces.  If  a shadow  shield 
geometry  is  chosen  for  applications  involving  human  or  very 
sensitive  equipment,  then  special  consideration  must  be 
given  to  the  regions  of  space  around  the  vehicle  that  might 
be  used  as  a working  area.  Figure  7-2  shows  the  four  main 
general  shield  configurations  that  might  be  selected.  As 
shown  in  this  figure  the  upper  limit  of  a shadow  shield  is 
the  four-pi  (4x)  shield.  This  design  provides  approximately 
the  same  amount  of  shielding  material  in  all  directions 
around  the  radiation  source.  In  addition  to  the  shadow 
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shield  and  four-pi  (4ir)  shield  configurations,  other  main 
configurations  include  the  two-pi  (2x)  shield  and  the 
preferential  (4x)  shield  shown  in  Figure  7-2. 

In  this  work,  NITAWL  is  used  to  generate  a 149  fine 
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Ultrahigh  Temperature  Vapor  Core 
These  preliminary  studies  are  performed  in 
order  to  find  shielding  designs  that  meet  performance 
specifications  for  acceptable  neutron  fluence  and  gamma  dose 
to  protect  the  sensitive  spacecraft  components  and  equipment 
(e.g.,  semiconductors,  tracking  systems,  etc.).  Fluence  and 
dose  scaling  or  normalization  for  the  burst  power  UTVR 
system  shielding  design  studies  is  achieved  by  assuming  a 
mission  time  of  one  hour  at  a core  power  level  of  1000  MW, 
(200  MWe) . The  UTVR  is  a highly  enriched  (>85%) , BeO 
externally-moderated,  circulating  fuel  reactor  with  UF,  as 
the  fissioning  fuel  (see  Chapter  2 for  more  information) . 
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geometry  with  S4  quadrature  using  the  XSDRNPM  code.  The 
materials  which  were  evaluated  include 
1.  natural  lithium  hydride  (LiH) , 


5.  (ZrHj  + 5%  B'°) , and 

6.  pure  tungsten  (W) . 

The  thicknesses  required  for  these  different  materials 


limit  of  10s  rad,  as  typical  radiation  damage  thresholds  for 
semiconductors,  are  shown  in  Table  (7-1) . All  of  the  11 

gamma  dose,  while  for  neutron  fluence  calculations  the  first 


neutrons,  LisH  has  the  largest  required  thickness  among 
the  Li  and  Zr  group  materials,  but  the  lowest  required  mass 
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per  unit  area.  Therefore,  Li6!!  is  the  best  examined 
material  for  neutron  attenuation.  By  adding  5%  B10  (by  atom 
weight)  to  the  zrH2  the  required  thickness  of  the  shield 
increased  slightly  (from  30.26  cm  to  31.10  cm)  but  the 
required  mass  per  unit  area  is  slightly  reduced  (from  1694.6 
to  1677.5  kg/mJ) . The  neutron  attenuation  ability  of  ZrHj 
configurations  on  a mass  basis  is  much  poorer  than  any  of 
the  LiH  configurations  (see  Table  7-1) . On  the  other  hand, 
tungsten  has  the  best  ability  to  attenuate  gammas  among  the 
examined  materials.  The  mass  per  unit  area  of  ZrHj 
configurations  to  reach  an  acceptable  level  of  both  gamma 
dose  and  neutron  fluence,  is  much  less  than  that  of  tungsten 
(because  tungsten  needs  a large  thickness  and  mass  compared 
with  the  other  materials  to  give  an  acceptable  level  of 
neutron  fluence) ; therefore,  ZrH2  is  the  most  efficient  of 
the  examined  materials  as  a space  reactor  shielding  material 
if  one  is  limited  to  a single  material.  Figures  7-3  to  7-7 
show  the  neutron  flux  behavior  as  a function  of  shield 
thickness  for  each  of  the  above  materials  except  tungsten. 
Each  of  these  figures  also  shows  the  minimum  thickness 
required  of  each  material  to  meet  the  acceptable  neutron 
fluence.  Figures  7-8  to  7-13,  on  the  other  hand,  show  the 
gamma  dose  versus  the  shielding  thickness  for  all  of  the 
examined  shield  materials.  Figures  7-11  to  7-13  also  show 
the  minimum  thicknesses  required  to  meet  acceptable  gamma 
dose  level  for  the  ZrHj,  (ZrHj+5%B'°)  and  tungsten  shields, 
respectively.  Tungsten  shows  the  best  ability  to  attenuate 


- Group  4 is  the  thermal  energy  group. 


Figure  7-3  Neutron  flux  versus  natural  LiH  shield 
thickness  for  1000  MWt  UTVR. 


- Group  4 is  the  thermal  energy  group. 

Figure  7-4  Neutron  flux  versus  Li7H  shield  thickness  for 


- Group  4 is  the  thermal  energy  group. 


Figure  7-5  Neutron  flux  versus  Li6H  shield  thickness 
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7-6  Neutron  flux  versus  ZrH 
for  1000  MWt  UTVR. 
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- Group  4 is  the  thermal  energy  group. 

Figure  7-7  Neutron  flux  versus  (ZrH,+5tB")  shield  thick- 
ness for  1000  MWt  UTVR. 


Figure  7-8  Gamma  dose  versus  natural  LiH  shield 
thickness  for  1000  MWt  UTVR. 
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Figure  7-11  Gamma  dose  versus  ZrH2  shield  thickness  for 


(pej) 


Figure  7-13  Gamma  dose  versus  pure  tungsten  shield 
thickness  for  1000  MWt  UTVR. 


gamma  rays  among  the  examined  materials  followed  by 
(ZrHj+5%B,0)and  then  by  ZrH2.  As  shown  in  Table  7-1,  the 
mass  per  unit  area  required  to  meet  acceptable  gamma  dose 
levels  is  1679.2  kg/m2  for  tungsten,  2077.6  kg/m2  for 
(ZrH2+5%B10) , and  2283.7  kg/m2  for  ZrH2.  On  the  other  hand, 
natural  LiH,  Li?H,  and  Li5!!  all  show  very  weak  gamma  ray 
attenuation  ability  due  to  their  low  Z values. 

Configurations  using  combinations  of  two  materials 
also  been  tested  to  determine  shield  thicknesses  required 
meet  acceptable  neutron  fluence  and  gamma  dose.  The  best 
material  from  the  single  material  shield  analysis  is  Li5!!. 
On  the  other  hand,  the  best  single  material  for  both  a 
neutron  and  gamma  shield  is  (ZrH2+5%B10) . These  two 
materials  are  each  used  with  tungsten,  the  best  single 
material  shield  for  gammas,  to  determine  optimum  shield 
configurations  that  will  minimize  the  shield  mass  and  stil 
meet  acceptable  levels  of  neutron  fluence  and  gamma  dose. 
Results  from  these  studies  are  shown  in  Figures  7-14  and  ~ 
15.  In  Figure  7-14  the  required  Li5!!  thickness  (i.e.,  the 
Li5!!  thickness  which  in  combination  with  the  tungsten 
thickness  gives  an  acceptable  neutron  fluence)  is  plotted 
versus  tungsten  thickness.  In  Figure  7-15  the  required 
(ZrH2+5%B’°)  thickness  (i.e.,  the 


thickness  of  ZrH2+5%B10 
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which  in  combination  with  the  tungsten  thickness  gives  an 
acceptable  neutron  £luence)  is  plotted  versus  tungsten 
thickness.  The  masses  per  unit  area  of  these  two-material 
shield  configurations  are  also  plotted  in  these  two  figures. 
On  Figure  7-14  the  required  thicknesses  for  the  combination 
of  U*H  and  tungsten,  that  give  an  acceptable  level  of  both 
neutron  fluence  and  gamma  dose  are  identified  by  point  A. 
These  thicknesses  are  found  to  be  8.27  cm  of  tungsten  and 
25.52  cm  of  Li6!!.  From  this  figure  the  shield  mass  per  unit 
area  of  this  33.8  cm  thick  combination  (Li*H  and  tungsten) 
is  found  to  be  1765.12  kg/m2.  From  Figure  7-15  the  required 
thicknesses  for  the  combination  of  (ZrHj+5%B'°)  and  tungsten 
that  give  an  acceptable  level  of  both  neutron  fleunce  and 
gamma  dose  are  also  identified  by  point  A.  These 
thicknesses  are  found  to  be  3.96  cm  of  tungsten  and  26.33  cm 
of  (ZrHj+5%B’°) . From  this  figure  the  shield  mass  per  unit 
area  of  this  30.3  cm  thick  combination  [ (ZrHj+5%B'°)  and 
tungsten]  is  found  to  be  2180.74  kg/m2. 

From  this  analysis,  ZrH2  is  the  most  efficient  of  the 
examined  materials  as  a space  reactor  shielding  material  if 
single  material  shield  configuration  have  been  used.  On  the 
other  hand,  the  combination  of  Li*H  and  tungsten  is  found  to 
be  the  most  efficient  of  the  examined  materials  as  a space 
reactor  shielding  material  for  a two  material  shield 
configuration . 
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CONCLUSION,  RECOMMENDATIONS, 
AREAS  OF  FURTHER  RESEARCH 


core,  externally  moderated  "cavity"  reactor  systems  fueled 
with  gaseous  UF4  using  the  S„-diffusion  theory  hybrid  model 


gas  density  distributions  in  the  core  are  investigated. 
These  two  cases  are  examined  to  determine  the  impact  of 


situation. 

A spherical  OCR  with  a 70  cm  core  radius,  and  a density 
factor  of  5.0  is  used  with  three  different  reflector 


thickness  on  the  performance  of  the  hybrid  model. 

Finally,  the  performance  of  the  hybrid  model  for 


Neutronic  calculations  are  performed  on  the  GCR  using 
the  hybrid  model  developed  in  this  work.  Different  GCR 
problems  are  examined  to  determine  the  performance  of  the 
hybrid  model  with  respect  to  accuracy  of  both  ke()  and  flux 
profile  predictions  and  acceleration  of  the  calculations. 
The  results  from  these  different  investigations  are  briefly 
summarized  below. 

performance  of  the  Hybrid  Model  for  Different 


density  factor 
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of  5.0  is  examined  (D.F.=5.0  corresponds 
to  the  fuel  atom  densities  given  in  Table  3-1) . It  is 
found  that  the  hybrid  model  overpredicts  kt)(  by  less 
than  1*  relative  to  the  standard  sn  approximation  when 
the  mathematical  interface  is  positioned  at  14  cm,  12 
cm,  and  10  cm  from  the  core-  reflector  interface  inside 
the  reflector.  The  hybrid  model  underpredicts  ke„  by 
less  than  1%  when  the  mathematical  interface  is 
positioned  at  8 cm  inside  the  reflector.  On  the  other 
hand,  it  is  found  that  the  hybrid  model  gives  large 
errors  when  the  mathematical  interface  is  located  at  6 
or  less  inside  the  reflector  for  this  particular 
conf iguration . 


The  reason  for  the  poor  performance  of  the  hybrid 
model  compared  to  that  of  the  standard  Sn  approximation 
when  the  mathematical  interface,  IMMP,  is  located  at  6 
cm  or  less  into  the  reflector  from  the  core-reflector 
interface  is  due  to  the  nonlinear  anisotropic  behavior 
of  the  flux  for  this  particular  density  factor  (i.e., 
D.F.=5.0)  at  these  positions. 

This  nonlinear  anisotropic  behavior  of  the  flux, 
which  continues  for  the  first  few  mean  free  paths 
inside  the  reflector,  renders  diffusion  theory  invalid 
for  the  same  distance  inside  the  reflector.  The  albedo 
values  employed  in  the  hybrid  model  are  generated  from 
the  converged  diffusion  theory  scalar  flux.  If 
diffusion  theory  is  not  valid  at  a certain  location,  it 


will  yield  incorrect  albedo  values  and  these  wrong 
albedo  values  will  cause  the  hybrid  model  to  give  poor 
results.  This  is  exactly  what  happens  when  IMMP  is 
located  at  6 cm  or  less  inside  the  reflector  for  this 
particular  configuration  (i.e.,  with  D.F.=5.0). 

It  is  found  also  that  the  hybrid  model  is  roughly 
five  times  faster  than  the  standard  Sn  approximation 
calculation  for  this  GCR  configuration, 
ii)  For  the  same  GCR  configuration,  the  core  flux 

behavior  predicted  by  the  hybrid  model,  when  the 
mathematical  interface  is  located  at  14  cm,  12  cm,  and 
10  cm  inside  the  reflector,  is  found  to  agree  exactly 
with  the  Sn  results;  the  core  flux  distribution 
predicted  by  diffusion  theory  differs  slightly  from  the 
hybrid  model  and  Sn  predictions. 

For  the  case  when  IMMP  is  located  at  8 cm  inside 
the  reflector,  the  results  from  the  hybrid  model  are 
still  good  but  start  to  deviate  slightly  from  that  of 
the  standard  sn  approximation. 

For  the  cases  when  IMMP  is  located  at  6 cm  or  less 
inside  the  reflector,  it  is  found  that  the  core  flux 
from  the  hybrid  model  is  not  in  good  agreement  with 
that  from  the  Sn  theory  approximation  for  this 
particular  configuration. 

For  mathematical  interface  locations  more  than  6 
cm  into  the  reflector  from  the  core-reflector 
interface,  it  is  found  that  the  hybrid  flux  (which  runs 


inside  the  reflector  renders  diffusion  theory  invalid  for 
the  same  distance  inside  the  reflector. 

The  albedo  values  used  in  the  hybrid  model  are 
generated  from  the  converged  diffusion  theory  scalar  flux. 

If  diffusion  theory  is  not  valid  at  a certain  location,  it 
will  give  wrong  albedo  values  and  these  wrong  albedo  values 
cause  the  hybrid  model  to  give  poor  results.  This  is 
precisely  what  happens  when  IMMP  is  located  at  6 cm  or  less 
inside  the  reflector  for  D.F.=5.0.  Also,  this  is  the  reason 
why  the  hybrid  model  gives  better  results  for  IMMP  located 
at  6 cm  inside  the  reflector  when  the  density  factor  is 
reduced  from  5.0  to  3.0. 

For  all  examined  density  factors,  the  hybrid  model  is 
found  to  be  about  five  times  faster  than  the  standard  S„ 
approximation . 

Pecf.PJ3iaag_e_s.f_the.Hybr.id.  Model  yitll.. Nonuniform 

Density  Distribution  in  the  Core. 

Two  different  idealized  and  extreme  cases  are 
investigated.  The  first  case  has  a density  factor  of  10.0 
in  the  central  region  of  the  core  and  a density  factor  of 
1.0  in  the  outer  region  of  the  core  while  for  the  second 
nonuniform  gas  density  distribution  case,  a density  factor 
of  1.0  is  used  in  the  central  region  of  the  core  and  a 

10.0  is  employed  in  the  outer  region  of 


density  factor  of 


the  hybrid 


ii) 


It  is  found  that  the  hybrid  model  is  roughly  five  times 


rgy  groups 


Conclusions 


i)  If  the  discrete  ordinates  approximation  is  applied  to 


The  developed  hybrid  model  uses  the  S„  approximation 


theory  is  used  for  the  rest  of  the  reflector. 

Therefore,  the  conclusion  that  can  be  drawn  from  the 
above  discussion  is  that  the  hybrid  model  advantage  should 
be  greater  for  GCR  systems  with  thicker  reflectors.  The 
presented  results  of  the  hybrid  calculations  performed  for 
GCR's  with  different  reflector  thicknesses  show  that  this  is 


core-reflector-shielding  configuration  calculations.  The 

model  over  the  Sn  method  is  probably  not  much  different  for 

reflector  system.  However,  since  the  total  computation  time 
required  by  the  coupled  gas  core-reflector-shield 
combination  can  be  quite  large,  the  total  savings  in  CPU 


ii)  The  developed  hybrid  model  uses  the  Sn  approximation  in 
the  gas  core  and  tor  the  first  few  mean  free  paths  in  the 
reflector.  Diffusion  theory  is  used  for  the  rest  of  the 
reflector. 

More  than  90%  of  the  hybrid  model  calculation  time  is 
used  in  the  Sn  approximation  section  for  the  gas  core 
region.  Since  it  has  been  shown  that  the  SK^  method  can  be 
significantly  faster  than  the  S„  approximation  (see  chapter 
4),  it  concluded  that  replacing  the  S„  approximation  section 
of  the  hybrid  model  with  the  SK,  method  could  significantly 
increase  the  advantage  of  the  hybrid  model. 

iii)  Space-time  dependent  calculations  are  needed  to 
determine  the  degree  of  coupling  which  can  exist  between  the 
neutron  field  and  the  gas  density  field  in  the  GCR.  Such 
space-time  calculations  performed  with  existing  time- 
dependent  Sn  transport  theory  codes  are  relatively 
expensive.  It  has  been  estimated  that  a one  to  two  order  of 
magnitude  speedup  in  computation  time  is  required  to  make 
such  calculations  practical.  It  is  concluded  that  the 
combination  of  the  hybrid  model  developed  in  this  work  with 
the  vectorization  and  multitasking  acceleration  techniques 
will  provide  over  a one  order  of  magnitude  improvement  in 
required  computation  time  (a  factor  of  15  to  25  is  estimated 
for  one-dimensional  problems) . It  is  also  concluded  that  a 
two  order  of  magnitude  improvement  could  be  approached  if 
the  Sn  approximation  is  replaced  by  the  SK,  method  in  the 


developed  hybrid  model. 


Areas  of  Further  Research 


combination  of  Hybrid  Model  With  Both  Vector jzation_an3 
Multitasking 

In  general,  it  is  found  that  the  hybrid  model  is  faster 
than  the  standard  S„  approximation  by  a factor  of  from  5 up 
to  8 times,  depending  on  the  problem  configuration  as 
explained  in  Chapter  6.  This  performance  of  the  hybrid 
model  may  be  further  improved  by  using  vectorization  (vector 
processing)  and/or  multitasking  (parallel  processing) . 

Vector  processing  is  a technique  for  reducing  the  amount  of 
processor  time  required  to  run  a program.  It  can  be  applied 
to  segments  of  FORTRAN  code,  such  as  a DO-LOOP  that  performs 
the  same  sequence  of  operations  on  successive  elements  of 
arrays.  Vector  processing  reduces  the  processor  time 
required  to  execute  such  loops  by  overlapping,  or  piplining, 
the  processing  for  the  array  elements.  Multitasking  is 
another  technique  for  improving  the  performance  of  computer- 
intensive fortran  programs.  Multitasking  can  be  used  on 
programs  that  can  be  divided  into  two  or  more  independent 
subtasks  and  run  in  a multiprocessor  configuration  (see 
Appendix  C for  more  information  about  vectorization  and 
multitasking) . 

If  the  code  in  which  the  hybrid  model  is  implemented 
(XSDRNPM-S  in  this  case)  is  written  and  arranged  in  a way  to 
take  advantage  of  vector  processing,  then  vectorization  will 
improve  the  performance  of  the  hybrid  model.  Parallel 
processing  can  also  be  of  advantage  to  the  hybrid  model. 


is  because  different  processors  can  be  used  at  the  same 
one  for  each  energy  group. 

The  combination  of  these  three  acceleration  techniques 
(i.e.,  the  hybrid  model,  vectorization,  and  multitasking) 
can  yield  a significant  reduction  in  CPU  computation  time 
over  that  of  the  standard  sn  approximation. 

Other  research  has  been  performed  on  discrete  ordinates 
using  vectorization  and  multitasking  techniques.  These 
investigations  show  that  the  acceleration  that  can  be 
achieved  by  using  the  vectorization  technique  is  roughly 
between  20%  to  30%  [87, 88);  the  acceleration  that  can  be 
achieved  by  using  the  multitasking  technique  is  about  2.5 
times  faster  than  the  standard  Sn  approximation  [89-91] . 

As  discussed  above,  the  hybrid  model  accelerates  the  Sn 
approximation  calculation  by  a factor  of  5 to  8 times 
faster.  If  parallization  accelerates  the  Sn  approximation 
calculation  by  a factor  of  2.5,  then  the  combined 
acceleration  factor  (for  both  the  hybrid  model  and 
parallization)  will  be  roughly  between  12  to  20,  depending 
on  the  GCR  configuration.  Moreover,  if  vectorization 
accelerates  the  Sn  approximation  calculation  by  20%  to  30%, 
as  explained  above,  then  the  total  combined  acceleration 
factor  (i.e.,  the  acceleration  of  the  hybrid  model  combined 
with  the  acceleration  of  both  vectorization  and 
multitasking)  will  be  roughly  between  15  to  25,  depending  on 


configuration. 


The  estimated  CPU  time  required  to  analyze  a 10-second 
GCR  transient  (space-time  dependent  dynamic  calculation) 
with  a CRAY  X-MP/48  supercomputer  is  more  then  1 hour  for  a 
four  group  one-dimension  Sn  calculation.  Therefore, 
according  to  the  combined  acceleration  discussed  above,  this 
CPU  computation  time  can  be  reduced  from  1 hour  to  roughly 
from  2 to  4 minutes.  This  is  a significant  reduction  in 
computation  time  for  these  calculations. 

The  performance  improvement  resulting  from  the  use  of 
these  three  techniques  should  be  even  more  significant  if 
used  with  two  dimensional  problems.  The  estimated  CPU  time 
required  to  analyze  a 10-second  GCR  transient  (space-time 
dependent  dynamic  calculation)  with  a CRAY  X-MP/48 
supercomputer  is  more  then  30  hours  for  a four-group  two- 
dimensional  S„  calculation.  Therefore,  the  reduction  in 
computation  time  will  be  significant  if  these  three 
techniques  are  combined  together  for  two-dimensional  space- 
time  dependent  dynamic  calculations. 

This  significant  reduction  in  computation  time  shows 
that  research  to  combine  these  three  acceleration  techniques 
is  very  important  and  recommended  as  an  area  for  further 
research. 

SK,-Piffusion  Hybrid  Model 

The  SK,  method  is  a high-order  transport  approximation 
introduced  for  solving  the  integral  transport  equation.  The 
method  relies  on  approximating  the  integral  transport 
kernels  by  a sum  of  diffusion-like  kernels.  The  integral 


equation  is  equivalent  to  a set  of  coupled  second-order 
differential  equations  for  which  blackbody  and  reflecting 
boundaries  are  established.  The  original  formulation  of  t 


SK,  approximation  has  a defect  when  applied  to  heterogeneous 
problems.  Therefore,  the  SK,  approximation  is  slightly 
altered  to  solve  the  integral  transport  equation  for 
heterogeneous  systems.  This  technique  can  also  be  applied 
to  problems  with  P,  scattering.  The  SK,  method  was  tested 
against  a variety  of  benchmark  problems  [81].  It  is  found 
that  the  SK,  method  is  up  to  16  times  faster  than  the  Sn 
approximation  for  some  cases.  It  is  noticed,  also,  from 
these  applications,  that  some  of  the  benchmark  problems  for 
which  this  method  proved  effective  have  properties  like 
>r  example,  in  one  of  the 


benchmark  problems  the  ratio  (£,/£,)  is  0.25.  This  ratio  is 
comparable  to  that  of  typical  GCR's.  Therefore,  it  is 
possible  that  the  SK,  method  might  be  superior  to  the  S„ 
approximation  for  GCR  calculations. 

An  SK,-dif fusion  hybrid  model,  analogous  to  the  S„- 
diffusion  theory  hybrid  model  developed  in  this  work  could 
significantly  reduce  the  CPU  time  required  for  GCR  neutronic 
calculations.  This  hybrid  model  would  depend  on  the  use  of 
an  initial  diffusion  theory  calculation  for  both  the  gas 
core  and  reflector.  Then  the  SK,  calculation  would  be  used 
for  the  gas  core  region.  The  steps  required  for  this  SK,- 
Diffusion  hybrid  model  are  similar  to  the  steps  employed  for 
the  Sn-diffusion  hybrid  model  presented  in  this 


dissertation.  Because  about  90%  of  the  total  computation 
time  in  the  developed  hybrid  model  is  spent  on  the  Sn 
portion  of  the  calculation  in  the  core  and  because  the  SK^ 
method  is  much  faster  than  the  S„  approximation  as  explained 
above,  an  SK,-diffusion  theory  hybrid  model  could  be 
significantly  faster  than  the  hybrid  model  developed  in  this 


Therefore,  development  o 


5 recommended  a: 


APPENDIX 


XSDRNPM-S:  A ONE-DIMENSIONAL  DISCRETE-ORDINATES 
CODE  FOR  TRANSPORT  ANALYSIS 

Introduction 

The  SCALE  version  of  XSDRNPM  which  is  called  XSDRNPM-S 
is  a one-dimensional  discrete  ordinates  transport  code  and 
is  the  latest  in  a series  of  codes  in  the  XSDRN  family.  As 
such,  it  contains  several  unique  characteristics,  though  a 
large  portion  of  the  theoretical  bases  and  intended  uses  of 
the  program  are  the  same  for  all  versions. 

XSDRNPM  is  an  acronym  for  2-gection  Eynamics  for 
Eeactor  nucleonics  with  Eetrie  Hodifications. 

XSDRNPM  is  a discrete  ordinates,  spectral  code  for  the 
generation  of  nuclear  multigroup  constants  in  the  fast,  and 
thermalization  energy  regions. 

The  code  is  written  to  accept  an  arbitrary  group 

inelastic  and  (n,2n)  scattering  matrices,  and  the  desired 
one-dimensional  reaction  arrays  such  as  absorption,  fission, 
(n,p),  etc.  The  present  Master  library  at  the  University  of 
Florida  contains  123  energy  group  data  extending  from  10'4 


Forward 


adjoint  solution  of  the  one-dimensional 


Boltzmann  transport  equation  is  performed  for  homogeneous, 
slab,  cylinder,  or  spherical  geometry.  This  solution  is 
performed  in  the: 

1.  multigroup  discrete  ordinates  approximation: 

2.  diffusion  approximation;  and 

3.  infinite  medium  approximation. 


Several  problem  ■ 


provided 


including: 


2.  k-calculation; 

3.  e-calculation  (flux  is  assumed  to  have  an  e*at  time 
variation) ; 

5.  outer  radius  search;  and 

6.  buckling  search. 

The  calculation  provides  spatial,  angular,  and  energy 
dependent  fluxes,  as  well  as  region  and  system  reaction  rate 
distributions. 

Spatially  and  energy  averaged  multigroup  parameters  are 
computed  for  subsequent  use  in  other  calculations.  Cell- 
averaged  as  well  as  region-averaged  microscopic  cross 
sections  are  available  for  each  nuclide  in  the  system.  A 
transport  cross  section  is  computed  and  edited  for  each 
nuclide  requested  from  the  library  tape.  The  multigroup 
parameters  can  be  output  on  tape  in  the  correct  format  for 
use  with  the  ANISN,  CITATION,  DOT,  EXTERMINATOR,  or  ROD 


programs,  and  the  output  routines  may 
provide  the  desired  format  for  other  c 
The  flexible  dimensioning  scheme 


easily  modified  to 


that  all  array  dimensions  are  set  for  the  particular  problem 
being  executed,  avoiding  the  wasted  core  storage  inherent  in 
normal  "fixed"  dimension  programs.  Furthermore,  if  all 
arrays  will  not  fit  in  core,  XSDRNPM  first  tries  to  gain 
additional  storage  by  putting  cross  section  data  on  an  I/O 
device;  if  there  is  still  insufficient  space,  any  fixed 
source  arrays  are  stored  externally,  and,  finally  the  flux 
moment  arrays  can  be  stored  externally,  if  need  be.  The 
complete  space  allocation  process  is  automatic  and  requires 
only  that  the  user  make  I/O  devices  available. 

XSDRNPM  differs  from  its  predecessor,  XSDRN,  in  several 
important  respects. 


* It  will  do  coupled  neutron-gamma  calculations. 

* It  allows  any  mixture  to  be  represented  to  an 
arbitrary  order  of  anisotropic  representation. 

* It  will  do  an  adjoint  calculation. 

* It  is  considerably  more  efficient  in  the  manner  in 
which  data  storage  is  used  and,  hence,  will  run  much 
larger  problems  in  less  core  storage. 

* The  resonance  calculation  is  removed  and  provided  in 
the  NITAHL  module  of  the  AMPX  or  SCALE  package. 

a It  employs  improved  thermal  flux  scaling  techniques 


problem  convergence. 


Input  specifications 


^-ordered , and 


defaults  have  been  provided  to  make  the  use  of  this 

* It  will  calculate  Sn  constants  for  any  even  order 
for  any  of  the  three  one-dimensional  geometries 
available. 

* Mixture-dependent  fission  spectra  are  calculated  and 

fissionable  nuclides  in  a problem. 

The  SCALE  version  of  XSDRNPM — XSDRNPM-S — has  several 
internal  differences  (mostly  transparent  to  the  user)  from 
the  AMPX  version,  the  major  one  being  the  use  of  double 
precision  arrays  for  fluxes  and  certain  other  key 
parameters.  The  use  of  double-precision  has  been  found  to 
circumvent  convergence  difficulties  encountered  in  some  ill- 
behaved  cases.  One  more  solution  option  for  calculation  is 
available  in  XSDRNPM-S  in  addition  to  the  three  above 
identified  options  of  XSDRNPM.  This  option  is  the  Bn  theory 
option  which  shares  many  of  the  same  restrictions  as  the 
infinite  medium  method.  However,  this  treatment  does  (as 
its  name  implies)  use  a buckling  approximation  to  account 
for  leakage  from  the  large  homogeneous  region,  thereby 
giving  higher  order  flux  moments  that  can  be  used,  for 
example,  to  determine  diffusion  coefficients. 


XSDRNPM-S  allows  a boundary  condition  to  be  specified 
for  each  of  the  two  "outside"  boundaries  of  its  one- 
dimensional geometries.  The  options  are: 

1.  Vacuum  boundary  - all  angular  fluxes  that  are 

directed  inward  at  the  boundary  are 


2.  Reflected  boundary  - 


boundary  of  slab,  T(»>0)=0, 
the  incoming  angular  flux  a 
boundary  is  set  equal  to  th 
outgoing  angular  flux  in  th 
reflected  direction;  e.g., 
left-hand  boundary  of  a sla 


3.  Periodic  boundary  - 


4.  White  boundary  - 


l|lJn(|l)  - t|f0„c(-ll)  . 
the  incoming  angular  flux  at  a 
boundary  is  set  equal  to  the 
outgoing  angular  flux  in  the  same 
angle  at  the  opposite  boundary, 
the  angular  fluxes  of  all  incoming 
angles  on  a boundary  are  set  equal 


boundary 


E «1M*. 
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. Albedo  boundary  - 


This  boundary  condition  is  used  as 
an  outer  boundary  condition  for 
cell  calculations  of  cylinders  and 
spheres  that  occur  in  lattice 


geometries. 

this  option  is  like  the  white 
albedo  boundary  condition  except 


dependent  albedo  multiplies  t 
incoming  angular  fluxes.  Thi 
option  is  rarely  used,  as  it 
difficult  to  relate  to  most 
practical  situations. 


Iteratisn .and  Convergence-Tests 

Two  parameters  are  used  to  specify  the  required  levels 
of  convergence  on  an  XSDRNPM-S  calculation.  These  are  EPS 
and  PTC,  both  given  in  the  5*  array.  The  flux  calculations 
proceed  through  a series  of  iterations  until  either 
convergence  is  achieved  or  the  specified  iteration  limit  is 


exceeded. 


The  basic  iteration  strategy  in  XSDRNPM  is  now 
described.  The  discrete  ordinates  difference  equation  is 
solved  for  the  first  angle  and  the  first  energy  group.  This 
sweep  generally  is  made  from  the  last  interval  boundary  to 

supplied  as  part  of  the  input  along  with  the  boundary 
conditions.  The  second  angle  is  then  calculated,  etc., 
until  all  angles  in  the  quadrature  are  treated.  At  the  end 
of  this  sweep,  new  scalar  fluxes  for  the  midpoints  of  all 
intervals  have  been  determined.  The  angular  sweep  continues 
until  either  the  point  scalar  fluxes  are  converged  to  within 
PTC  or  until  the  code  makes  IIM  inner  iterations.  An 
exception  to  this  "inner  iteration"  pattern  occurs  on  the 
first  outer  iteration  (defined  below)  whenever  a fission 
density  guess  is  used,  instead  of  the  flux  guess.  In  this 
case,  the  program  uses  one-dimensional  diffusion  theory  to 
determine  a scalar  flux  value  for  all  intervals  and  the 
angular  sweeps  are  not  made  until  the  second  outer 
iteration.  After  the  first  group  is  completed,  the 
calculation  goes  to  the  second  group  and  repeats  the  above 
procedure.  This  continues  until  all  groups  have  been 

The  pass  through  all  groups,  angles,  and  intervals  is 
called  an  outer  iteration.  Host  of  the  convergence  checks 


involving  all  energy  groups  and  are  made  against  the  EPS 
parameter  mentioned  above.  For  a coupled  neutron-gamma 
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iteration,  the  following  checks 


4.  |1.0  - **1  & EPS, 

5 . R 1 1 . 0 - | i EPS,  and 

6.  B|1.0  - | * EPS. 

R is  a convergence  relaxation  factor  and  is  set  internally 
to  0.5  in  XSDRNPM-S.  If  all  convergence  criteria  are  net, 
if  I CM,  the  outer  iteration  maximum,  is  reached,  or  if  ITMX, 
the  maximum  execution  tine,  is  exceeded,  the  problem  will  be 
terminated  with  full  output!  otherwise,  another  outer 
iteration  will  be  started. 


Cross  Section  Weighting 

The  collapsed  cross  sections  from  XSDRNPM  are  written 
as  an  AMPX  weighted  library  on  logical  unit  3 which  can  be 
used  directly  by  the  nodule  to  do  transport  calculations. 
They  can  also  be  output  in  cards  image  or  in  a binary  format 
suitable  for  the  ANISN,DOT,  or  MORSE  code.  Four  weighting 
options  are  provided: 

Cross  sections  are  generated  in  a manner  consistent 
with  mocking  up  a cellular  configuration  as  a homogenized 
region;  VIZ.,  the  spatial  "disadvantage  factors"  are  taken 


section  weighting. 


sections  are  defined  by  the  following  equation: 

W of"  bgV-  dz  ml)  f^dE$(E,  f)  o (E.  z)  (A-9) 


N(r)  is  the  number  density  of  the  element  at  r, 


0 is  the  cell  average  flux  in  group  g defined  by: 


f dE  f dz  4>(E,f) 
d>  . lS!2 I22U 

/„„* 


(A-10) 


where  V is  the  cell  volume,  a (E , r)  is  the  element  cross 
section  for  the  process  at  E and  r,  0(E,r)  is  the  neutron 
flux  per  unit  energy  and  volume  at  E and  r,  N is  the 
homogenized  number  density  defined  by: 

f dZ  W<?) 

N - , (A-ll) 

dz 


L, 


and  oflce"  is  the  cell-averaged  cross  section  . 

The  fine  group  cross  sections  input  into  XSDRNPM  are 

averaging  codes  assume  some  characteristic  flux  shape,  such 


as  1/E  or  fission  spectrum, 


the  effects  o 


p structure  is  generally  chosen 
averaging  spectrum  are  not 


important.  The  fluxes  calculated  by  XSDRNPM  are  integrated 
(in  energy)  values.  Therefore,  the  internal 
definition  of  a cell-averaged  cross  section  is: 

- g Nj  £ LVt  £ *,(J) o„U)  <»-“> 


N,  is  the  number  density  in  zone  j, 
i is  a spatial  interval  index, 

4V,  is  the  volume  of  interval  i, 
g is  a fine  group  index, 

os(j)is  the  group  g value  of  the  microscopic 
IH  is  the  total  number  of  intervals.  Also, 


4 >,(i)  - dE((>(B,f)  , (A-14) 
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s produced  that  is  weighted 
over  each  material  region  in  which  a nuclide  occurs.  This 
set  of  cross  sections  is  given  according  to  the  following 
equation: 


Oa""1  <t£  Vj  - g bVi  Y,  $g<.i)og{j)  (a- 18) 

the  volume  of  zone  j.  This  scheme  produces  a set  of  cross 
sections  for  a nuclide  for  every  zone  in  which  it  occurs, 
iii.  ^Region"  weighting 

These  cross  sections  are  averaged  only  over  the  regions 
which  contain  the  material.  Therefore,  this  kind  of 
weighting  produces  one  set  of  cross  sections  for  a nuclide 
but  weighted  over  a composite  spectrum  made  up  of  all 
spectra  from  zones  where  the  nuclide  is  present. 


dz  Mii)  Jao  dflfr  (E,  i)  a (E,  i) 
f dz  N{i)  |ao  dEo(E.i) 
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"region"  averaging  scheme  produces  one  set  of  cross 
r nuclide. 

4.  "inner"  tell  weighting 

regions  selected  by  the  user  in  the  input.  This  allows  one 
to  mock-up  a calculation  which  can  include  non-zero  leakage 
at  the  outer  boundary  of  the  cell;  e.g.,  an  explicit  cell 
mock-up  can  be  described  within  a homogeneous  description  of 


Transport  cross  section 

XSDRNPM  will  produce  transport  cross  sections  according 
to  two  different  models.  The  difference  between  the  two 
lies  in  the  method  used  to  provide  spatial  weighting. 

A point  definition  of  the  transport  cross  section  may  be 


defined 


(A— 21) 
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diffusion  coefficient, 


perpendicular  buckling, <Pt  i 


Definition  2 

Average  fluxes. 


e used  as  weighting  factors. 
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The  terms  involving  the  integral  on  the  right  side  of 
Equation  (A-21)  are  multiplied  by  "weighting"  factors  a 
specified  by  Definition  1 or  Definition  2.  Define  the 
second  term  on  the  RHS  of  Equation  (A-21)  as 


vil. 


This  yields  the  following  equations  f 
Definition  1 


h definition. 
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The  current  (J)  in  Equation  (A-26)  and  Equation  (A-27) 
is  approximated  by  the  leakage  spectrum  defined  in  Equation 
(A-22)  and  the  fine-group  transport  cross  section  is 
calculated  by  subtracting  Equation  (A-26)  from  Equation 
(A-23)  or  Equation  (A-27)  from  Equation  (A-24),  depending  on 
the  definition. 

The  fine-group  transport  cross  sections  are  then 
averaged  over  a leakage  spectrum,  I»9  , calculated  as 

Lg  - DgBg  V * dE  f dA  f*'  dO  0 •<(>  (E,  I.  Q)  {A_28, 

The  first  term  of  the  RHS  in  Equation  (A-28)  represents 
the  perpendicular  leakage,  while  the  seoond  term  represents 
the  parallel  leakage.  In  the  above  equation,  dA  is  a unit 
surface  area  on  the  outside  of  the  system,  and  $(E,r,n)  is 
the  flux  per  unit  energy,  volume,  and  solid  angle  at  E,  r, 


(A-29) 


The  definitions  previously  given  for  a transport  cross 
section  are  modified  for  "zone"  or  "region"  weighting 
consistent  with  the  schemes  discussed  in  the  "cross  section 
weighting"  section. 


APPENDIX 


XSDRNPM-S  SUBROUTINE  DESCRIPTION  AND 
BRANCHING  STRUCTURE 


(B-l) , while  Figure 


Main  Flow 


The  main  Program  (named  MAIN)  calls  SETUP  to  read  the 
first  data  block  and  indirectly  calls  XSDRIV  through  ALOCAT. 

several  10  buffers  and  reads  the  first  block  of  data.  XSDRIV 


1-1  DRTRAN  - DRTRAN  is  the  control  routine  for  most  of  the 


ilizatic 


SET  UP 


r nixing, 


d storage  required  f 
ck  the  information  read  in  the  first 
ock.  Subroutine  MIX  is  then  called  t 
e cross  sections.  Subroutine  PLSNT  i 


called  to  read  the  rest  of  the  input  data  and 
EDIT  is  called  to  edit  the  group  boundaries 


MIX  controls  the  cross  section  mixing.  It 
prints  the  mixing  table  and  initializes  each 
mixture  to  zero.  It  then  reads  the  input 
cross  section  library  and  mixes  the  nuclides 
according  to  the  mixing  table.  It  calls  JLLl 
to  generate  a pointer  vector  for  the  2D  cross 
section  arrays  and  calls  MIXEM  to  do  the 
actual  mixing.  After  all  the  mixing  is  done, 
it  loops  back  through  the  mixtures,  calling 
DIFCOF  to  initialize  the  diffusion 
coefficients  for  each  mixture.  It  calls  NCHI  to 
normalize  the  mixture  fission  spectrum  and  to 
finish  generating  the  mixture  diffusion 
coefficients.  During  this  loop  the  mixture 
cross  sections  are  printed  if  requested  by  the 
input.  Finally  the  variable  array  pointing  to 
the  beginning  of  each  mixture  array  is 
calculated. 

MIXEM  mixes  the  cross  sections  from  the  input 
library  into  a mixture.  If  it  is  called  for 


one-dimensional  cross  sections,  it  mixes  them 
into  the  front  of  the  P0  cross  section  array. 
If  the  problem  is  an  adjoint  problem,  the  cross 
sections  are  mixed  in  adjoint  order. 

DIFCOF  - DIFCOF  initializes  the  mixture  diffusion 

coefficients  to  the  within-group  P0  scattering 
cross  section  for  the  mixture. 

NCHI  - NCHI  normalizes  the  mixed  fission  spectrum  for 
a mixture,  generates  a transport  cross  section, 
and  a mixture  diffusion  coefficient. 

PLSNT  - PLSNT  reads  the  third,  fourth,  and  fifth  data 
blocks  and  sets  up  the  storage  pointers  for  the 
rest  of  the  problem.  If  enough  storage  is  not 
available  for  the  problem,  data  are 
automatically  scheduled  to  be  held  on  random 
access  until  the  problem  either  will  fit  or  no 
more  data  are  available  to  be  stored  on  random 

STORXS  - storxs  loads  the  mixture  cross  section  into 
memory  or  onto  random  access  as  required. 

ADJNT  - ADJNT  is  used  to  invert  the  sources,  the  flux 
guess,  the  white  boundary  albedos, and  the  broad 
group  definition  array  for  an  adjoint  problem. 

STAPE  - STAPE  writes  the  input  fixed  sources  to  random 


33#  Array 


storage 


onto  random  access 
the  fluxes  are  not  held  in  memory. 

1-11  REDF  - REDF  reads  a flux  guess  from  a previously 
written  scalar  flux  tape. 

1-12  FILL  - FILL  is  a utility  routine  used  to  fill  an  array 
with  a constant  supplied  in  the  calling 
sequence. 

1-13  DOQ  - DOQ  computes  a default  quadrature  set  for  the 
problem. 

1-14  ROOTS  - ROOTS  solves  for  the  roots  of  the  orthogonal 
polynomials  when  generating  the  discrete 
ordinates  quadrature  set. 

1-15  FIND  - FIND  locates  the  zeroes  of  a polynomial  by  a 
range-halving  technique.  It  is  used  in  the 
default  quadrature  generation  process. 

1-16  QUADWT  - QUADWT  evaluates  the  weights  corresponding  to 
the  roots  of  the  appropriate  Legendre  polynomial 
for  use  in  generating  the  discrete  ordinates 
quadrature  set. 

1-17  PNCNST  - PNCNST  checks  the  discrete  ordinates  quadrature 
set  and  computes  the  Pn  constants  that  convert 
between  the  angular  flux  and  the  flux  moments. 

1-18  GETSTS  - GETSTS  computes  the  number  of  sets  of  weighted 
cross  sections  to  produce  for  a given  nuclide, 
and  computes  an  atom  density  for  each  set. 


DFLTBN 


banding  arrays 


the  input  data. 


from  the  group  boundaries 


rse  velocities 
print  out  the 


energy  and  lethargy  group  bounds.  It  then 
prints  out  the  mixture  by  zone  array,  the  order 
of  scattering  by  zone  array,  the  activity 
material  and  reaction  type  arrays,  and  the 
quadrature  weights,  direction  cosines,  reflected 
direction,  and  weight  times  cosine  arrays. 
Finally  it  prints  out  the  constants  for 
converting  between  angles  and  moments. 


Branch  2 


2-1  OUTERS  - OUTERS  is  the  flow  control  routine  for  the  flux 
and  search  calculation.  It  calls  routines  to 
generate  geometry  dependent  factors  needed  in 
the  discrete  ordinates  flux  calculation,  to 
normalize  the  fixed  source,  if  necessary,  to 
calculate  and  normalize  the  fission  source,  to 
calculate  the  scattering  source,  to  calculate 
the  flux  according  to  discrete  ordinates, 
diffusion  theory,  infinite  homogeneous  medium, 
or  Bn  theory,  and  to  check  convergence  and 
compute  the  search  parameter.  OUTERS 


controls  the  loops  over  energy  groups  that 
define  an  outer  iteration  loop.  These  are 
actually  two  nested  loops,  one  over  bands  and 
another  over  groups  within  a band. 

The  outer  iteration  monitor  table  is  printed 
at  the  end  of  each  outer  iteration.  This  table 
shows  the  current  value  of  the  eigenvalue  and 
the  search  parameter  as  well  as  the  status  of 

ZMSH  - ZMSH  finds  the  first  and  last  spatial  intervals 

GEOMF  - GEOMF  computes  the  geometry  dependent  factor 
needed  in  the  discrete  ordinates  flux 
calculation.  It  checks  that  the  input  spatial 
mesh  is  nonnegative  and  in  increasing  order.  It 
modifies  the  mesh  for  a search  problem. 

FIXSRC  - FIXSRC  computes  the  total  fixed  source  per  group 
and  normalizes  the  source  to  the  input 
normalization  factor. 

RETRVE  - RETRVE  retrieves  a block  of  the  fixed  source 
data  from  random  access  when  needed. 

REPLCE  - REPLCE  writes  the  fixed  source  onto  a random 
access  unit  when  this  is  needed  because  of 
insufficient  memory  or  because  it  was  requested 
by  the  input. 

FISSRC  - FISSRC  computes  the  fission  source  before  each 


integral 


integral  of  the  previous  source.  Finally, 
fissrc  scales  the  fluxes,  flux  moments,  boundary 
fluxes,  and  the  fission  densities  by  X. 

BASADR  - This  subroutine  is  used  to  compute  the  base 

address  array  (IBA)  pointing  to  the  beginning  of 
each  P,  piece  of  each  mixture  when  the  cross 
sections  are  stored  on  random  access  storage 
rather  than  in  memory. 

POSTN  - POSTN  computes  the  pointers  to  the  next  group  of 
cross  sections. 

SCATSC  - SCATSC  computes  the  scattering  source  for  a 

group.  It  also  calculates  82  needed  for  flux 
calculations  performed  by  Bn  theory. 

BN  - BN  calculates  the  fluxes  and  moments  according 
to  Bn  theory.  B2  for  the  calculation  is  computed 
by  subroutine  SCATSC.  BN  also  estimates  a 
leakage  and  returns  it. 

TRISOX,  - TRISOL  is  used  to  solve  the  tridiagonol  matrix 
in  Bn  theory  for  the  fluxes. 

CELL  - This  subroutine  computes  fluxes  using  an 
infinite  homogeneous  medium  approximation. 

DT  - DT  performs  the  diffusion  theory  flux 

calculation  required  when  a fission  density 
guess  is  supplied,  or  when  flagged  in  the  input. 

BDYFLX  - This  subroutine  is  used  to  initialize  the 
boundary  fluxes  after  the  fluxes  have  been 


initialized  by  a diffusion  calculation  using 
fission  density  guess. 


2-16  INNER  - INNER  is  the  routine  that  does  the  discrete 

ordinates  angular  flux  calculation.  It  uses  the 
current  scalar  flux  and  flux  moments  to 
calculate  the  scattering  source  for  the  current 
group.  It  then  loops  over  the  angles, 
calculating  the  angular  flux  at  each  mesh 
boundary  using  a weighted  diamond  difference 
scheme,  traversing  the  mesh  in  the  direction  of 
the  current  angle.  During  the  sweep,  it 
computes  updated  scalar  flux  moments.  When  it 
reaches  a system  boundary,  it  saves  the  new 
boundary  flux. 

2-17  REBALN  - REBALN  generates  flux  rebalance  factors  for  an 
inner  iteration. 

2-18  XCEL  - XCEL  is  used  to  accumulate  data  used  in  the 

outer  iteration  flux  acceleration  calculation. 

2-19  SCALFC  - SCALFC  calculates  the  matrix  used  in  upscatter 
scaling  and  then  calls  UPS CAL  to  solve  for  the 
scaling  factors. 

2-20  UPSCAL  - UPSCAL  solves  the  upscatter  scaling  matrix  for 
the  scaling  factors. 

2-21  EXCEL  - FXCEL  does  outer  iteration  flux  extrapolation 
when  this  option  has  been  triggered  in  the  2$ 
Array.  Much  of  the  programming  in  FXCEL  is  ad 


hoc  with  little  testing,  and  this  option  is  not 
used  by  any  of  the  present  SCALE  procedures. 
convrg  - CONVRG  is  called  at  the  end  of  every  outer 

iteration  to  determine  if  the  code  has  converged 


search  problem. 


Branch  3 

output  calling  Seouence 

OUTPUT  - OUTPUT  is  where  the  scalar  flux  is  printed, 

punched,  and/or  written  to  an  output  file  when 
requested.  The  broad  group  fluxes  are  also 
printed  and  saved  here  if  requested. 

INVERT  - INVERT  is  used  to  flip  fluxes  and  flux  moments 
in  an  adjoint  problem  back  to  forward  order  for 
printing  and  writing  on  the  flux  tape.  If  the 
fluxes  are  output  on  random  access  storage,  they 
are  also  rewritten  there  in  forward  order  for 
later  use  in  cross  section  weighting, 
activities,  and  balance  tables. 

PRTF  - PRTF  is  used  to  print  the  scalar  fluxes. 

WRTF  - WRTF  is  used  to  write  the  scalar  flux  tape. 

BFLUX  - BFLUX  collapses  the  fine  group  fluxes  to  broad 
groups  and  writes  them  on  the  end  of  the  flux 


BT  galling  Sequence 


BT  - BT  computes  the  pointers  for  collecting  the 

balance  tables,  clears  the  arrays  that  will  be 
used,  writes  a head  record  if  the  balance  table 
information  is  being  saved,  and  calls  SUMARY  to 
compute  and  print  the  balance  tables. 

SUKARY  - SUMARY  calculates  the  information  printed  in  the 
balance  tables  from  the  angular  and  scalar 
fluxes  and  the  cross  sections. 

DWOT  - DWOT  is  a double-precision  version  of  subroutine 
WOT  used  to  print  out  the  angular  fluxes. 

BASADR  - This  subroutine  is  used  to  compute  the  base 

address  array  (XBA)  pointing  to  the  beginning  of 
each  p(  piece  of  each  mixture  when  the  cross 
sections  are  stored  on  random  access  storage 
rather  than  in  memory. 

RETRVE  - RETRVE  retrieves  a block  of  fixed  source  data 


n needed. 

N computes  the  pointers  to  the 


TOTAL  - TOTAL  is  used  to  accummulate  the  totals  over 
groups  and  over  zones  in  the  balance  tables. 
PRTBAL  - PRTBAL  is  used  to  print  the  balance  tables. 
WRTBAL  - WRTBAL  is  used  to  write  the  balance  table 
information  to  an  external  file. 


- SUMFEW  collapses  the  balance  tables  to  the  few 
group  structure  for  printing  if  requested. 

Branch  5 

ftCTY  Calling  Sequence 

This  subroutine  computes  the  pointers  needed  for 
computing  the  reaction  rates,  clears  the  arrays 
which  will  be  used,  and  calls  subroutine  RRATE 
to  compute  the  reaction  rates. 

RRATE  calculates  the  reaction  rates  specified  in 


FEWG  controls  the  cross  section  weighting 
calculation.  It  sets  up  the  pointers  for  the 
calculation,  calls  RECRD1  to  output  the  title 
record  on  the  weighted  cross  section  tape,  sets 
up  the  random  access  storage  to  be  used  in  the 
weighting,  call  CURNTS  to  calculate  the 
weighting  functions  to  be  used,  calls  RECRD2  to 
output  the  group  structure  records  on  the 
weighted  tape,  calls  WATE  to  do  the  cross 
section  weighting,  and  calls  CFLUX  to  print  the 
cell  averaged  fluxes,  the  flux  disadvantage 
factors,  and  cell  averaged  currents. 


RECRD1  - RECRD1  computes  the  values  needed  tor  the  first 
record  of  the  weighted  cross  section  library  and 
then  writes  the  record. 

GETSTS  - GETSTS  computes  the  number  of  sets  of  weighted 
cross  sections  to  produce  for  a given  nuclide, 
and  computes  an  atom  density  for  each  set. 

CURNTS  - CURNTS  computes  the  zone  fluxes  and  currents 
that  will  be  used  as  weighting  functions  in 
producing  a weighted  cross  section  library. 

BASADR  - This  subroutine  is  used  to  compute  the  base 

address  array  (IBA)  pointing  to  the  beginning  of 
each  P,  piece  of  each  mixture  when  the  cross 
sections  are  stored  on  random  access  storage 
rather  than  in  memory. 

POSTN  - POSTN  computes  the  pointers  to  the  next  group  of 

RECRD2  - RECRD2  writes  the  group  boundary  records  onto 
the  weighted  cross  section  library.  It  computes 
an  average  inverse  velocity  and  fission  spectrum 
for  the  problem  and  writes  these  at  the  end  of 
the  neutron  group  boundary  record.  It  then 

WATE  - WATE  does  the  cross  section  weighting 
calculations. 

WT1D  - WT1D  does  the  broad  group  weighting  of  the  one- 
dimensional reaction  cross  sections. 


TRNPRT  - TRNPRT  applies  the  PI  tern  to  the  transport 
cross  section  in  the  weighting  calculation. 
CFLUX  - CFLUX  computes  and  prints  the  flux  disadvantage 
factors  as  well  as  the  zone  fluxes  and  currents. 


Branch  7 

SPOUT  calling  Sequence 

SPOUT  - SPOUT  outputs  the  library  format  asked  for  after 
a weighting  calculation  has  been  done. 

INFACE  - INFACE  controls  the  flow  of  converting  a 

weighted  cross  section  tape  to  CCCC  interface 
format.  The  variables  for  the  first  two  records 
on  the  interface  file  are  set  up  and  the  records 
are  written.  The  nuclides  on  the  weighted 
library  are  then  looped  over  and  those  selected 
for  the  interface  file  are  put  in  the  correct 
format  and  written  on  a scratch  file  while  the 
information  for  the  head  record  for  the  nuclide 
is  accumulated.  The  head  record  is  then 
written  on  the  interface  file  and  the  rest  of 

to  the  interface  file. 

TR1D  - TR1D  loads  the  one-dimensional  cross  section 
weighted  library. 


format  library 


COLAPS  - COLAPS  is  used  when  a CCCC  ISOTXS  cross  section 
interface  file  is  requested  to  collect  the 
correct  one-dimensional  cross  sections  to  be 
placed  on  the  interface  file. 

TRSCT  - TRSCT  is  used  to  load  the  transfer  arrays  for 
the  CCCC  format  library  from  the  weighted 


APPENDIX 


ysstprUa.t.ian  and  Multitasking 

Vectorization  and/or  Multitasking  are  two  techniques 
which  nay  be  used  to  accelerate  FORTRAN  codes  [87-90] . 

Vector  processing  is  a technique  for  reducing  the 
amount  of  processor  time  required  to  run  a program.  It  can 
be  applied  to  segments  of  a FORTRAN  code,  such  as  a DO-LOOP 
that  performs  the  same  sequence  of  operations  on  successive 
elements  of  arrays.  Vector  processing  reduces  the  processor 
time  required  to  execute  such  loops  by  overlaping,  or 
piplining,  the  processing  for  the  array  elements. 

The  VS  FORTRAN  Version  2 compiler  can  automatically 
generate  vector  instructions  in  appropriate  places  in  the 
object  module  as  it  compiles  the  VS  FORTRAN  (language  level 
66  or  77)  source  program  (this  process  is  called 
Vectorization) . However,  this  compiler  examines  each  D0- 
LOOP  and  selects  for  vectorization  only  those  statements 
that  can  be  safely  and  effectively  transformed.  Statements 
that  cannot  be  safely  vectorized,  or  would  not  run  faster  in 
vector  mode,  are  compiled  into  standard  scalar  instructions. 

The  compiler  selects  certain  of  the  analyzed  loops  as 
being  eligible  for  vectorization.  This  selection  is  based 
on  cost  estimates  (in  CPU  cycles)  of  executing  each  loop  in 
vector  and  scalar  mode. 


compile-tine  option  VECTOR  invokes 


vectorization  process,  which  produces  programs  that  use  the 
IBM  3090  vector  Facility.  The  facilities  in  the  VS  FORTRAN 
Version  2 compiler  used  for  vectorization  include: 

i.  coapile-llme  option  and  Suboptione  f°r  vectorisatjpn: 

The  compile-time  option  VECTOR  invokes  the 
vectorization  process,  which  produces  programs  that  use  the 
IBM  3090  Vector  Facility.  The  following  suboptions  are 
provided  by  the  VS  FORTRAN  Version  2 Compiler: 
a.  LEVEL  (0/1/2) 

a-1  LEVEL  (0)  specifies  no  vectorization. 
a-2  LEVEL  (1)  vectorizes  entire  DO-loops, 

provided  all  statements  within  are 
eligible.  If  a portion  of  a loop 
is  not  eligible,  none  of  it  is 
vectorized.  This  option  yields 
faster  compilation  and  less  code 
rearrangment  than  LEVEL  (2) , but  it 


may  produce  less  vectorization. 

) vectorizes  individual,  eligible 
statements  in  a DO-loop,  even  if 


the  entire  loop  is  not 
vectorizable.  As  much  of  a loop  as 
possible  is  translated  into  vector 
object  code  and  the  remainder  is 
translated  into  scalar  code.  This 
produces  the  best  execution-time 


performance,  but  does  take  more 
compilation  time  and  does  result  in 
some  code  rearrangement. 

REPORT  ( LIST/XLIST/TERM) 

This  option  produces  a vector  analysis 
report.  The  report  lists  the  results  of 
the  compiler's  source  code  analysis  and 
shows  which  DO-loops  have  been  vectorized 
and  which  will  execute  in  scalar  mode. 

The  report  can  be  used  to  examine  the 
source  program  for  possible  coding 
refinements.  Certain  sections  that  were 
not  vectorized  may  be  recoded  to  make  them 
eligible  for  vectorization. 
b-l  LIST  produces  the  simple  format  of  the 
vector  report  in  the  output  listing 


XLXST  produces  the 
vector  report 
file. 


extended  format  of  the 
in  the  output  listing 


terminal  only. 
INTRINSIC/NOIHTRINSIC 


Specifying  intrinsic,  authorizes  the 
compiler  to  select  vector  intrinsic 
functions  in  vectorizable  statements. 


REDUCTION/  NONREDUCTION 


Specifying  reduction  permits  vectorization 
of  reduction  operations. 

The  SIZE  suboption  permits  the  user  to 
instruct  the  compiler  to  generate  vector 
code  for  a specific  section  size. 

2.  Intrinsic  Functions  (Vector  Mathematical  Functions) : 
Vector  versions  of  most  of  the  VS  FORTRAN  Version  2 
intrinsic  mathematical  functions  are  provided  with  the 
product.  Thus,  calculations  using  these  functions  within 
DO-loops  are  eligible  to  take  advantage  of  the  speed  of 
vector  processing.  Because  these  vector  mathematical 
functions  have  the  same  names  as  their  nonvector 

compiler  automatically  selects  the  correct  version. 

The  basic  unit  of  vectorization  by  the  VS  FORTRAN 
Version  2 Compiler  is  the  DO-loop.  The  compiler  analyzes 
only  DO-loops  for  vectorization.  For  example,  loops  formed 
by  a GO  TO  statement  to  a target  preceding  it  will  not  be 
vectorized.  The  loop  must  be  formed  by  a DO  statement  to  be 
vectorized.  The  range  of  a DO-loop  consists  of  all  of  the 
executable  statements  that  appear  following  the  DO  statement 
that  specifies  the  DO-loop,  up  to  and  including  the  terminal 
statement  of  the  DO-loop.  If  a DO  statement  appears  within 
the  range  of  a DO-loop,  we  have  a nest  of  DO-loops.  In  this 
case,  the  compiler  analyzes  the  inermost  DO-loop  in  a nest 


of  DO-loops,  and  selects  the  single  loop  whose  vectorization 
will  lead  to  the  fastest  execution  of  the  entire  nest.  As 
with  many  other  vectorization  actions,  there  are  ways  to 
write  the  statements  in  the  nest  to  increase  the  compiler's 
chances  of  exploiting  possible  vectorization  opportunities. 
The  compiler  also  make  an  economic  analysis  before  it  start 
vectorization.  Consider,  for  example,  the  following  nest  of 


X(X,J)  = AA(I, J)  + BB(I,J) 

15  Y(I,J)  = AA(X,J)  * BB(I,J) 

The  compiler  considers  both  DO-loops  as  vectorization 
candidates.  An  economic  analysis  is  made  with  the  available 
information  on  the  loop  limits,  dimensions,  stride,  cost  of 
instructions,  and  so  forth.  Depending  on  the  economic 
analysis,  the  compiler  might  vectorize  the  outer  loop,  the 
I-loop  here,  rather  than  the  inner  loop  on  J.  The  compiler 
will  also  investigate  the  possibility  that  the  nest  will 
execute  fastest  if  neither  loop  is  vectorized.  It  is  quite 
possible  that  scalar  execution  might  be  faster  than  vector 
execution.  If  this  is  the  case,  the  compiler  will  generate 

Any  job  that  uses  the  IBM  3090  vector  facility  needs  to 
be  run  on  a computer  that  has  the  necessary  hardware  to 
perform  vectorization.  Also,  at  the  Northeast  Regional  Data 
Center  (NERDC) , it  is  required  that  a SETUP  VECTOR  statement 


be  included  in  the  job  JCL  statement.  The  vectorization 
SETUP  and  EXEC  statements  are  illustrated  below. 


/•SETUP  VECTOR 

//  EXEC  F0RTVCLE,SUBLIB1='SYS1.ESWLIB', 

//  OPTIONS5® ' VECTOR ( LEVEL ( 2 ) REPORT  SIZE(LOCAL) ) MAP  LC(80)  • 


Multitasking  is  another  technique  for  improving  the 
performance  of  FORTRAN  programs.  Multitasking  is  designed 
to  improve  turnaround  time  for  computer-intensive 
applications.  The  amount  of  turnaround  improvement  that  can 
be  achieved  depends  on  the  application.  The  primary  factors 
affecting  turnaround  performance  are  the  number  of 
processors  and  the  percentage  of  the  application  processing 
that  can  be  performed  in  parallel.  It  can  be  used  on 
programs  that  can  be  divided  into  two  or  n 


subtasks  and  run  on  a multiprocessor  configuration.  This 
allows  the  subtasks  to  be  executed  concurrently,  when  the 
system  workload  permits,  thereby  potentially  reducing  the 
turnaround  time  required  by  the  job. 

The  VS  FORTRAN  Multitasking  Facility  (MTF)  is  a feature 
of  the  VS  FORTRAN  library  that  takes  advantage  of  the 
multitasking  capabilities  of  the  MVS/Extended  Architecture 
(MVS/XA)  and  MVS  operating  systems.  MTF  may  be  invoked 
using  several  user-accessible  subroutines  supplied  by  the  VS 
FORTRAN  Version  2 library  and  a subparameter  on  the  EXEC 


statement. 


MTF  provides  the  following  functions: 

1.  Initialize  the  MTF  environment 

2.  Schedule  parallel  subroutines  for  execution 

3.  Synchronize  completion  of  the  parallel  subroutines. 

The  VS  FORTRAN  Version  2 Library  creates  the  MTF 

environment  required  to  run  the  separate  parts  of  a program 
in  parallel.  The  initialization  occurs  when  the  keyword 
AUTOTASK  is  specified  in  the  PARM  parameter  of  the  EXEC 
statement  used  to  run  the  program.  The  AUTOTASK  keyword 
specifies  two  subparameters.  The  first  subparameter  is  the 
name  of  a load  module  that  contains  the  program's  parallel 
subroutines.  The  second  subparameter  is  the  number  of 
subtasks  that  should  be  created  for  the  program. 

One  of  the  MTF  subroutines,  NTASKS,  returns  the  number 
of  subtasks  specified  with  the  AUTOTASK  keyword  in  the  PARM 
parameter  of  the  EXEC  statement  for  the  job  step.  NTASKS 
may  only  be  called  from  a main  task  program. 

The  main  task  program  schedules  a parallel  subroutine 
for  execution  by  calling  the  MTF  subroutine  DSPTCH.  The 
parallel  subroutine  is  assigned  to  a subtask  and  may  be 
executed  in  parallel  with  the  main  task  program,  or  other 
instances  of  itself,  or  other  parallel  subroutines.  The 
main  task  program  can  schedule  multiple  instances  of  a 
parallel  subroutine  for  parallel  execution  by  repeating  the 
call  to  DSPTCH  using  the  same  parallel  subroutine  name,  but 
passing  different  arguments  with  each  call. 


The  main  task  program  synchronizes  its  own  execution 
with  the  completion  of  the  parallel  subroutines  by  calling 
the  MTF  subroutine  SYNCRO.  SYNCRO  causes  the  main  task 
program  to  wait  until  all  of  the  currently-scheduled 
parallel  subroutines  complete  their  execution. 

A parallel  subroutine  is  coded  as  a FORTRAN  subroutine 
that  follows  certain  rules  required  for  proper  operation 
with  MTF.  In  addition  to  data  independence,  the  rules  are: 

1.  A parrallel  subroutine  cannot  call  the  MTF 
subroutines  NTASKS,  DSPTCH,  and  SYNCRO.  Such 
calls  can  only  be  used  in  the  main  task  program. 

2.  A parallel  subroutine  must  not  use  common  blocks. 
All  communication  with  the  main  task  program  must 
be  by  way  of  data  that  is  received  through  the 
subroutine's  dummy  argument  list. 

3.  The  only  FORTRAN  I/O  statements  allowed  in  a 
parallel  subroutine  are  PRINT  statements  or  WRITE 
statements  that  are  directed  to  the  error  unit. 


This  unit  is  usually  unit  6,  but  may  have  been 
changed  when  the  vs  FORTRAN  Version  2 library  w 


n alternate 


installed. 

The  dummy  argument  1 
cannot  contain  a 
(asterisk) . 

The  parallel  subroutine  may 
series  of  subroutines  that  c 


a parallel  subroutine 
return  specifier 

actually  be  coded  as  a 
sail  one  another.  All 


these  subroutines  operate 


subroutine's  subtask  environment  and  must  follow 
the  rules  of  a parallel  subroutine.  However,  a 
subroutine  that  is  called  within  this  environment 
can  use  alternate  return  specifiers. 

A parallel  subroutine  is  always  invoked  in  its 
last-used  state.  If,  for  example,  a parallel 
subroutine  has  initialized  a variable  with  a DATA 
statement,  then  that  variable  has  that  value  the 
first  time  that  copy  of  the  parallel  subroutine  is 
used.  Should  that  value  be  modified,  then  that 
modification  is  available  to  the  next  execution  of 
that  copy  of  the  parallel  subroutine.  It  cannot, 
however,  control  which  copy  of  a parallel 
subroutine  is  used  when  the  parallel  subroutine  is 
scheduled.  Therefore,  this  parallel  subroutine 
must  not  depend  upon  residual  values  from  previous 


APPENDIX 


GENERAL  DESCRIPTION  AND  FLOW  DIAGRAM  OF 
OUTERS  SUBROUTINE 


In  this  appendix,  the  flow  diagram  of  the  OUTERS 
subroutine  is  presented  in  Figure  D-l.  OUTERS  is  the  flow 


to  the  user  desired  approximation.  For  example,  it  calls 
the  DT  subroutine  to  calculate  the  flux  according  to  the 


definitions  for  most  of  the  OUTERS  subroutine  variables  that 


given  ir 


Flow  diagram  of  the  OUTERS  subroutine. 


(Continued) 


Figure  D-X  (Continued) 


Figure  D-l  (Continued) 
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Figure  D-X  (Continued) 


Figure  D-l  (Continued) 
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Table  D-l  (continued) 


Pointwies  convergence 


Accumulated  n 


APPENDIX 


GENERAL  DESCRIPTION  AND  FLOW  DIAGRAM 
OF  THE  DT  SUBROUTINE 


diffusion  theory  flux  calculation  required  when  a fission 
density  guess  is  supplied,  or  when  flagged  in  the  input. 
This  flow  diagram  shows  the  main  steps  performed  when  the 

listing  of  the  DT  subroutine  is  also  shown  in  this  appendix. 
Finally  definitions  for  most  of  the  DT  subroutine  variables 


Figure  E-l  Flow  diagram  of  the  DT  subroutine. 
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(Continued) 


Figure  E-l 
(Continued) 


(Continued! 


(continued) 


APPENDIX 


GENERAL  DESCRIPTION  AND  FLOW  DIAGRAM  OF 
the  INNER  SUBROUTINE 

The  INNER  subroutine  performs  the  discrete  ordinates 
angular  flux  calculation.  It  uses  the  current  scalar  flux 
and  flux  moments  to  calculate  the  scattering  source  for  the 
current  group.  It  then  loops  over  the  angles,  calculating 
the  angular  flux  at  each  mesh  boundary  using  a weighted 
diamond  difference  scheme,  traversing  the  mesh  in  the 
direction  of  the  current  angle.  During  the  sweep,  it 
computes  updated  scalar  flux  moments.  When  it  reaches  a 
system  boundary,  it  saves  the  new  boundary  flux.  The  flow 
diagram  of  the  INNER  subroutine  is  shown  in  Figure  F-l. 

The  Fortran  listing  for  this  subroutine  is  also  shown  in 
this  appendix.  Finally  the  definitions  for  most  of  the 
variables  which  appear  Figure  F-l  are  given  in  Table  F-l. 


Initialization 
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(continued) 


n space  boundary  f 


n angle  boundary  £ 


5 -SECTION 


ENERGY  GROUP  STRUCTURE  OF  CROSS 

LIBRARIES  USED  IN  THIS  WORK 

this  Appendix  the  energy  group  boundaries  1 
nd  4 neutron  cross-section  libraries  are 
, G-2,  and  G-3  respectively. 


Group  energy  boundaries  (ev) 

From 

i . 105E+7 

1.000E+7 

8 . 187E+6 

7.408E+6 

5.488E+6 

12 

14 

15 

16 

3.012E+6 

17 

3.012E+6 

20 

2 231E+6 

2 019E+6 

1 . 653E+6 

23 

1 496E+6 

1 353E+6 

1 225E+6 

1.225E+6 

1.108E+6 

ible  G-l  (Cc 


286 


p boundaries  (ev) 


[1) 


[«) 


[7] 


«S?S 


■ATSS1. 


E.  C.  bittern  and  B.  ^ 

TSSZfik'  <£2”i;5)."“ Band  corporatio" 

iasKlr 


[31] 


over  System  with  a Flssioning-Plasma  Disk 
* Proc.  ot  14th  IEEE  International 
lasma  Science,  Washington,  DC  (Jun. 


■r;“£s;r 


mmmmm 


ss™~“ 


[91] 


BIOGRAPHICAL  SKETCH 


Mohammed  K.  Alfakhar  was  born  on  July  1,  1952,  In  Baghdad 
Iraq.  In  1972,  Mohammed  Alfakhar  matriculated  to 
Almustansirai  University  in  Baghdad  and  received  the  degree  of 
Bachelor  of  Science  in  Physics  in  1977.  He  continued  his 
graduate  study  in  Nuclear  Physics  at  the  University  of  Baghdad 
in  Baghdad  and  received  the  Master  of  Science  degree  in  1979. 

Mohammed  Alfakhar  spent  the  following  seven  years  doing 
research  and  teaching  at  the  University  of  Baghdad.  In 
Spring,  1988,  Mohammed  Alfakhar  entered  the  graduate  program 
in  nuclear  engineering  at  the  University  of  Florida  in 
Gainesville,  Florida. 

Mohammed  Alfakhar  is  married  to  the  former  Rwaida  Alsidi. 
Their  daughters  are  Zenah,  13  years  old  and  Sura,  6 years  old. 
Their  sons  are  Ahmmed,  10  years  old  and  Kassim,  4 years 


